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SECTION 1 


INTRODUCTION 

1.1 Definition of the Problem 

The evaluation of the aerodynamic pressure is an important 
tool for the design of aeronautical and space vehicles. Current 
methods do not satisfy the requirements of generality, flexi- 
bility and efficiency.* A general theory of potential aero- 
dynamic flow around a lifting body having arbitrary shape and 
motion is presented here. The theory is based upon the classical 
Green theorem approach, and provides a tool for the evaluation 
of the aerodynamic pressure acting on the surface of the body. 
Comparison with existing results shows that the proposed method 
is not only more general and flexible, but also at least as fast 
and accurate as existing ones. The concept of the Green function, 
is fundamental to the theory and therefore an ample discussion 
of this concept is given in Section 2, where the classical ex- 
pressions for the subsonic and supersonic Green function are 
derived using a novel approach which, it is hoped, gives a clear 
physical interpretation for these expressions. A detailed out- 
line of this report is given in Subsection 1.3. In the following, 
a short analysis of the method (lifting surface theory) currently 
used in the design of aircraft is presented. 


* 


An excellent analysis of the state of the art is given by 


25 


Ashley and Rodden. 




The problem of the evaluation of the pressure acting on 

a surface of a body immersed in a fluid stream has always 

attracted the attention of the scientists. In particular, 

interest in the theory of unsteady potential aerodynamics, 

which is a basic tool in dynamic aeroelasticity , has been 

12 3 

frowing steadily in the last fifty years. ' ' Until re- 
cently, the attention of the researchers has been concen- 
trated on lifting-surface theories in which the body is 
assumed to have zero thickness; the solution of the problem 
is reduced to an integral equation relating pressure and 
downwash or similar quantities. An excellent analysis of 
the recent literature in this field is given in Refs. 2 
and 3 . 

Two major objections can be raised about the lifting- 
surface theories. First, the numerical solution of the 
problem is rather complicated. The difficulties are re- 
lated to the complicated form of the kernel function and 
the numerical integration of improper integrals. An ex- 
cellent analysis of the numerical problems which are en- 
countered in the various lifting-surface formulations is 
given in Ref. 3. Attempts to circumvent these difficulties 

4 

have been presented more recently, but the present situation 
can be considered still unsatisfactory. 


la 



The second objection is that the lifting-surface 
theory cannot be easily generalized to include more com- 
plicated geometries and motions. Geometries which need 

further exploration are for instance the effect of thickness 

25 

of the wings, and the wing-body interference. The results 
shown in Ref. 5 show that the lifting-surface theory is 
being used beyond the limit of its validity and that, in 
some cases, the results are still unsatisfactory. Similar 
is the situation for motions which do not fall into the 
categories of either harmonic oscillation or impulsive start. 

An attempt to circumvent the present "impasse" situation 
is described in Ref. 6, where a three-dimensional body of 
arbitrary geometry which executes an arbitrary motion in an 
incompressible fluid is considered. The method has the 
advantage of satisfying the boundary conditions on the true 
location of the surface of the body. This is a considerable 
improvement with respect to the lifting surface theories 
since it does not involve the use of zero thickness. Thus, 
complicated geometries can be easily analyzed. Furthermore, 
the formulation is suitable for treating arbitrary motions. 

On the other hand, the formulation still involves a singular 
kernel, which implies the above-mentioned difficulty of 
evaluating the integrals in the principal sense. In addition, 
the method is limited to incompressible flow, since the 
Green theorem for the Laplace equation is used. 


2 



1.2. Formulation of the Problem 


In the following, the isentropic inviscid flow of a 
perfect gas, initially irrotational , is considered. Under 
this hypothesis, the flow can be described by the velocity 
potential 4* . The equation of the unsteady aerodynamic 
potential, given by Garrick^ - is 

** v 2 <± . , , 

^ Dt l 1,1 

where a is the speed of sound, V 1 is the Laplacian operator, 
and 

= - 2 _ + Vi-v 1.2 

T>t 

is the total time derivative (the subscript "c" reminds 
that V<f> should be treated as a constant in order to ob- 
tain the second total time derivative) . Consider a frame 
of reference such that the undisturbed flow has velocity 
u. in the direction of the positive x-axis. Then the 
speed of sound is given by 



Furthermore, it is convenient to introduce the perturbation 
potential , such that 

<f> * U - ( x4 f) 1.4 

Note that 0 in the undisturbed flow. Combining Eqs . 

1,1, 1.3 and 1,4 yields the equation for the perturbation 
potential 


3 



1.5 


< v 2 <p = U 1 [6* <f.)‘ f„ +■ f/ ft, 

+ 2 U 1 <f,, +(U<p x )<p, f f r 

t 2 U, [(Uf.)<f x . * f,f,. 

+ (r ') i [f t + u <p, ♦ i f; + f,*)] v 2 f 


For the sake of simplicity, it is convenient to separate 
the linear terms from the nonlinear ones; thus Eq. 1.5 
is rewritten as 



1.6 


where 



A. z. A + u 0 

cU 04 


1.7 


is the linearized total time derivative and the nonlinear 
terms are given by 





4 



This is the equation for the unsteady potential compressible 
(subsonic or supersonic) aerodynamic flow. 

A very general approach is considered here by assuming 
that the body immersed in this flow has arbitrary shape 
and is moving with arbitrary motion. Thus, the surface of 
the body is represented in the general form 

5 3 (x,f l, t)= O 1.9 


where the subscript B stands for Body. The boundary con- 
dition on the body is given by 


= O 1.10 

X) t 


or 

+ vj> ■ = O 1.11 

By using Eq. 1.4, Eq . 1.11 yields 



^ Wt 





1.12 


or 

0$ = 1.13 

with 




JUt 


1. 14 


Furthermore, as mentioned above, the boundary condition at 
infinity is given by 



1.15 


Finally, the pressure can be evaluated from the Bernoulli 


5 



Theorem 


r 



or the linearized Bernoulli Theorem 



which yields, for the pressure coefficient 


c _ p - f>~ _ _ 2 / J. i.is 

p ' j o u l \ v. at »x / <W 

rw «* 


Note that no assumption was made on the Mach number 

M » Uj 1.19 

Thus, the above equations are valid for both subsonic and 
supersonic flow. 

1.3. Method of Solution 

The method of solution presented here is based upon 
the well known Green function technique. The Green func- 
tions for the linear unsteady subsonic and supersonic flow 
are derived in Section 2, and used in Section 3 to derive 
an integral representation of the potential cf> at any point 
in the field (control point) in terms of the values of 
and /fin on a surface surrounding the body and the wake. 
In Section 4, it is shown that the integral representation 


6 



can be simplified under the assumption of small perturbation; 
and even further simplified for small vibration around a 
configuration fixed with respect to the frame of reference. 
Next, in Section 5, it is shown how lifting surface theories 
can be obtained as limiting cases of the formulation pre- 
sented here. 

Finally, in Section 6, the problem of small vibration 
around a fixed configuration in subsonic flow is analyzed in 
detail. It is shown (Appendix C), that if the control point 
is on the surface of the body, the problem reduces to an in- 
tegral equation. The question of existence and uniqueness 
is also discussed in Section 6 and, in particular, it is 
shown that, for a zero- thickness flat wing, the integral 
equation operator becomes singular. 

Hence, in order to verify the limit of applicability 
for low values of the thiclcness ratio, the steady subsonic 
flow around very thin wings is solved numerically. The re- 
sults are presented in Section 7, whereas the conclusions 
are discussed in Section 8. 
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SECTION 2 


THE GREEN FUNCTION FOR SUBSONIC AND SUPERSONIC 
LINEAR UNSTEADY POTENTIAL FLOW 

2.1. Introduction 

The Green function for the linear unsteady potential 
flow for the whole space (which represents a unit-impulsive- 
source) is the solution of the problem 


- 1 2.i 

or oJl 

06 

with 

Cj = o dt 2,2 

In Eq. 2.1, S is the well known Dirac delta function de- 
fined by 

04 ao 

JJff&Afc " JF/ F (* vlf * 1 5 2 - fc - 1 -) 2 3 


The solution of Eq. 2.1 for the subsonic case (subsonic 
Green's function) is given by 


<$ _ - —L J (fc t - fc+T) 

where is the usual Dirac delta function and 

. i 

r p = I (■‘-‘Ap* 

and 


T , 




tAU- k,) 

- 


2.6 


8 



where 


f= ^i-n j 2 - 7 

8 9 

Although this result is well known, ' for the sake 
of completeness, Eq. 2.4 is derived here (Subsection 2.3) 
using a particularly instructive procedure. 

9 10 

Similarly, the supersonic Green function is given by ' 



I 

4-zjr^ 



2.8 


with 



2.9 


and 



M (*-*, ) ± 


r 

a. 


2.10 


where 

3 - J 2.ii 

Equation 2.8 is derived in Subsection 2.4. 

2.2. Galilean Transformation 

In order to obtain the Green function of the linear 
unsteady potential flow, it is convenient to use a Galilean 
transformation, such that the new frame of reference is 
rigidly connected to the undisturbed flow. Then the differ- 
ential equation reduces to the well known wave equation for 
which the Green function is well known. Then, by using 


9 



the inverse transformation, the Green function for the un- 


steady linear potential flow is derived. 

In order to. avoid transformation of generalized func- 
tions (distributions) it is convenient to consider the non- 
homogeneous equation of the linear unsteady potential * 




2.12 


This equation reduces to the nonhomogeneous wave equation 
for a coordinate system g rj T rigidly connected to the 
undisturbed flow. The Galilean transformation relating 
the two systems is given by 


*--5* u .. T 1 = 1 2=5 t ‘ z 


2.13 


whose inverse is given by 


gz x- u t 

Since 

rj = l| 

2.14 

9 _ ^ ^ *0 

/ 3_ / 0 -0 _ © 

2.15 

’ 0x 9r)' 

95, 9* = 9t “3* 


Equation 2.12, in the 

new frame of reference, reduces to 


<*£ + ^ * 2 ^ 

- ’ = T e ($ t «S,z\ 

2.16 





with 


r .(4, 1.5. = Q (% 


O T 

oi i 






2.17 


* 

Note that Q is a fictitious prescribed source distribution, 
while F in Eq. 1.6 depends on the unknown itself. 
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The solution of Eq. 2.16 is given by 


11 


p/*- 




2.18 


where 



( 5 - 5 .)*] 

- J 




and 



T (L , 1> ?. > T .j 


with 

r, • z - p /«. 

or 


2.19 


2.20 


2.21 



2.22 


Equation 2.18 is equivalent to saying that the Green function 
for the wave equation in the space is"^ 



2.23 


Finally, in order to obtain the solution of Eq. 2.12, it is 
sufficient to express Eq. 2.18 in terms of the original 
variables x, y, z and t. It should be noted however, that 
the inverse transformation is not given by Eq. 2.14, since, 
according to Eq. 2.21, one has 


11 



2.24 


T-T. * P /*- 

Thus, the inverse of the transformation 

*-*. = + (*-T. ) - + ^f 2.25a 

H- 4 . VI. 2 ‘ 25b 

Z - 2, = 2.25c 

t-t, = X -T, ^ p/*„ 2 . 25d 


has to be considered. The inverse is given by (*) 


5-5, -- 1 , i 

1 - M 2 - 

^X-X, + Mr) 

2.26a 

V 1, °M. 


2.26b 

-- *-Z, 


2.26c 

T-T, = p /^ = 

1 1 [ ? * M(x-X,)] 

«« 1-M 1 L J 

2 . 26d 


(*) Equation 2.25a yields 

[u-o-a-o] 2 - = «* [(*,-?.)* ♦(i-i.Ms-s.r] 

By solving for £,-^, , one obtains 

= — (*-*,) t Hr 
I - 

that is, Ecr. 2.26a. Combining this equation with Eq. 2.19, 
yields, in agreement with Eq. 2.26d 


f 


0 -mV 
. j [M (*-*,) i V 2 " 


(*-*,)* i 2Mr(x-*,) r 1 


^V^VCz-a ,) 1 
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with 




0-^)[(h-h,? 



2.27 


The interpretation of the double sign which appears in Eq. 
2.26a is discussed later in this section. For clarity, it is 
convenient to consider the subsonic case and the supersonic 
one independently (Subsections 2.3 and 2.4, respectively). 
Note that, in either case* 


A 

r - 


M(Vi) + p 


2.28 


and that the Jacobian of the transformation is given by 


j. >H. i 3.) _ 

3(!., >).,«.) 


which shows that 


unSi' i Jk Mli 

P P P 




O 




2.29 


2.30 


By combining Eqs . 2.25 and 2.27 yields 
r z ^ (x-* ( )' + ( i-M*) Of-'O'+C *-*) 2 

= { \-l,Y + 2 + * Y 

s ($-§,)* + 2 Mp ($-•£.)+ n*U-U 4 + 

* [tf - p ] 2 


13 



2.3. .Subsonic Green's Function 


Consider the subsonic case. 

and 

I*-*.) <7 


For M < 1 , one has r = 

2.31 


Thus, Eq. 2.26d, by eliminating the absolute-value sign, 
may be rewritten as 



2.32 


By substituting this equation and Eq. 2.26a into Eq. 2.25a, 
one obtains 


x-x, 



p[± «(*-*,) + r^J 


2.33 


which is satisfied only if the lower sign is used. 

Thus, by using the lower sign, the inverse trans- 
formation for the subsonic case is given by 


with 


This 

2.18 


S-l. = p M ->] 

1 - 1 , - H-'i' 

ZL-l, -- 


2.34 


T - x ' - P A - 5 ^p[ r p- M(k 'H 

. p t n (4-4,) 2-35 

transformation may be used in order to express Eq. 
in the frame of reference (x, y, z, t) . This yields 


i'M) m IT 
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2.36 


= ///— <2/ 

JJI 4/7p I 6. : 6 


-?K 




where 


with 


fJI^ f Q J T dx '^ dz ‘ 

[°J T = 


2.37 



2.38 


Equation 2.36 shows that the Green function for subsonic 
case is 


<5 


— — S(t,-t + T) 

4/i*> v ' 


P 


as shown in Eq . 2.4. 


2.39 


2.4. The Supersonic Green Function 

Consider the supersonic case. For M > 1, Eq. 2.25a 

yields 

X-X, >0 2.40 


which shows the well known property of supersonic flow 
that any point can influence only the points downstream. 
Equation 2.40 implies 

(*-*< ) > r B 2.41 

Thus, Eq. 2.26d may be rewritten as 


] 


2.42 
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2 2 

Note that B = M - 1 > 0. Finally, combining Eq. 2.25a, 
and 2.42 yields 


X-X. 




2.43 


which is satisfied by both signs. 

The interpretation of this double sign lies in the 
well known fact that a disturbance A, (Fig. 1) traveling 
at supersonic speed, influences a given point B two times. 

The first by "backward-traveling waves"coming from position 
A + , and the second by "forward- traveling waves" coming from 
position A . In other words, the supersonic inverse trans- 
formation is not a one-to-one but a one-to-two transformation. 
Thus, for supersonic case, the inverse transformation is 
given by 


with 


1 - 1 , 

s-c, 




+ Mr, 


2 - 2 , 


t-*. * f7«- = [*(*-*■) i r.J 

o4 


2.44 


r _ 

E> 



2.45 


This transformation may now be used in order to express 
Eq. 2.18 in the frame of reference (x, y, z, t) . This 
yields, in analogy to the subsonic case. 
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2.46. 


lf(x, 3 . 2.0 



I 

4nr s 



‘M*/, A 


where 


wi th 



2.47 



2.48 


Equation 2.46 shows that the supersonic Green function is 
given by 



_ I 

4 »f, 


■5 (t,-t * . £(t,-t+T_) 


in agreement to Eq. 2.8. 


2.49 
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SECTION 3 


GENERALIZED HUYGENS' PRINCIPLE FOR THE NONLINEAR EQUATION OF 
THE UNSTEADY AERODYNAMIC POTENTIAL 

3.1. Introduction 

As mentioned in the introduction, the purpose of this 
analysis is to obtain a representation of the potential in 
terms of its values (and the values of its derivatives) on the 
surface of the body and the wake. A representation of this 
type for the wave equation is called Huygens' Principle 1 ^ 

(or Kirchhoff's formula). The corresponding formula for the 
unsteady compressible (subsonic or supersonic) aerodynamic 
potential will be called Generalized Huygens' Principle . 

In order to obtain this principle, it is convenient to 
follow a procedure similar to the one used in Ref. 11 to de- 
rive the Huygens' Principle. There are two major differences 
between the two problems. The first is that Eq. 1.6 is more 
complicated than the wave equation examined in Ref. 11. The 
second is that the surface is assumed to be changing in time. 
This assumption is necessary if enough generality is desired, 
in order to include arbitrary motions, like the roll for in- 
stance. In order to simplify the generalization of the pro- 
cedure used in Ref. 11, it is convenient to make use of the 

. 7 

theory of distributions developed by Schwartz. 

In order to do this, it is convenient to introduce a few 
basic definitions. 
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Note first that the equation of the aerodynamic potential 
given by Eq. 1.6, is not valid on the wake, where discon- 
tinuities on c|> exist. Thus, consider the volume V in which 
Eq. 1.6 is valid. At any instant of time, this volume is 
given by the whole physical space except the volume, V fi , 
occupied by the body and an infinitesimally thin layer, V w , 
representing the wake. Define the function E (see Fig. 2) 


E (x ,y , z , t) = 1 on V 

= 0 otherwise 


3.1 


This function represents the domain of validity of the 
equation of the potential and will be called "domain function”. 
Consider the surface of discontinuity of the function E, that 
is the surface, Zl , surrounding the volume . Let 


S(x, Z, t)= 0 3.2 

be the equation of the surface X 

Note that the surface 2 is composed of two branches. 

The first, is the surface of the body given by Eq. 1.9. 

The second is the surface, 2.^, of the wake 

= ° 3,3 

Note that this surface Z. , is considered twice, since Z is 
a closed surface. In other words, the upper side and the lower 
side of the wake are considered to be two independent surfaces 
having the same equation (but opposite outwardly-directed 
normal) . 
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Finally, for later convenience, the four-dimensional 
gradient of the function E is introduced. Consider first the 
four -dimensional outwardly-directed* normal to the surface 
y 1 , defined by 


- OS 
v = 

\nSl 


( 3S/2x ) 
1 ) ds/dy 

la s { ] 3S/3* ( 
l -95)/3t ; 


3.4 



3.5 


It should be noted that, in writing Eq. 3.4, it has been assumed 

implicitly that the vector D S is equidirected with the four 

— * 

dimensional normal V . This condition can always be satisfied 
(by a suitable change of sign in Eq. 1 if necessary). 

Next, note that along the direction of the normal V , 
the function E behaves as a step-function. Thus its directional 
derivative in the direction of the normal v is a Dirac delta 
function on the surface X , which will be indicated with the 
symbol 



It may be worth noting that the integrals of <5^ are equivalent 
to hypersurface integrals (the surface H is on the four di- 
mensional space). 


* "Outwardly" is understood as "going from the body into the 
fluid", that is, from the region E = 0 into the region E = 1 
(see Fig. 1) 
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3.7 


JJJf F ^ -- <jfj> Fdl 


By using Eq. 3.6, the four-dimensional gradient of the 
function E can be written immediately, since Eq. 3.6 is equi- 
valent to 


OE = iiJ-. S T JLL 

| nsj 


3.8 


where Eq. 3.4 has also been used. 

By separating the spatial components from the time com- 
ponent, Eq. 3.8 yields 

VE VS _!_ 

21 ias| 

3.9 

d £ _ J 35 _j_ 

db |0$J 

Note that 


4 !L -- Qel* u = J L 

dt £ dt |0S| 


3.2. Green's Theorem for the Equation of Aerodynamic Potential 
In order to obtain the generalized Huygens' principle, it 
is convenient to follow the general method that leads to the 
Green theorem. Multiply the equation of the aerodynamic 
potential (in the form given by Eq. 1.6) by the Green function 
G and subtract Eq. 2.1 (definition of the Green formula) 


multiplied by 


f ' 


« frf - i ^ - £ % )= 5 F- 

where the arguments of cp and its derivatives are x^, y^, 


3.11 


z.^ and t 1 , while the arguments of G and its derivatives are 


x i ■ x ' ^i ~ y- 


- z, and t^ - t. 
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Making use of the identities 


y. (aVb)= Va-tf.fc + 3.12 


and 


d ( rx clt> \ - da. djo + a d l l> 

d£ t ' dJtJ dt, £tfc, cU, 1 3.13 


Equation 3.11 reduces to 


^.(qy(p-fys)-J. £(g = SF-f<T 


3.14 


Multiplying Eq. 3.14 by the domain function E, defined by 
Eq. 3.1, and integrating over the whole four-dimensional 
space yields 


§ £ [ y ■(* - r~ 6 )- f )] 

■4 

= jlll E (GF -tfSjM.M, 


3.15 


where the subscript 1 indicates the dummy variable of in- 
tegration , and 


l . £ t ui 

dJt, fr, “ 3*. 


By suitable integrations by parts, Eq. 3.15 yields 

■4 

-JUJ J, ^ 

. {ft! £ (GF-fS)dVM, 

- -4 
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or, by using Eqs . 3.9 and 3.10 


•M[ v ,S-(G9,f-<fVfi} 



= JJ1I E (SF -t pi)dV, di, 

- mS 



3.18 


which is the Green theorem for the equation of the aerodynamic 
(subsonic and supersonic) potential. Note the presence of the 
factor which shows that the integral on the left hand 

side of Eq. 3.17 is a surface integral, as it is indicated by 


Eq. 3.7. 

Finally, making use of the definition of the Dirac delta 
function (Eq. 2.3) yields 

Ecp = jjjj £6 F 


if [v.s Vfi ) -± g * -&|4 I nsf-'^dt, 3 . 
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3.3. Generalized Subsonic Huygens' Principle 

In this subsection, Eq. 3.19 is specialized to the subsonic 
case for which the Green function is given by Eq. 2.4, which 
may be rewritten as 



with 

a T . 

where T is given by Eq. 2.38. Note that 


3.20 


3.21 



and 


_d| T 

<stt 


, 0*. / 7 ^ tl V 9x, / 


3.23 


Combining Eqs . 3.19 to 3.23 yields 

4n £ cp = - JJIf £ 1 S f F rfv/.rtt, 

rv;s^«P - - — ^ 1 _L <5_ J I ns| «w,<& 

l ' <**, <**. J Tp * T 1 

- 06 


3.24 


A, 


Note that the integrands of the integrals on the right hand 
side contain products of distributions (either <5^ S r or 
<* t W T /0t, ). For the sake of convenience, integrals of this 
type are discussed in Appendix A where it is shown that, for 
any "regular" function f, 


|os t^dt, . f[ ( ] T i^r 


3.25 


and 


*h f -£ |[f] T |vsr^ T 


dt 


Z T 


3.26 


T . 


In Eqs. 3.25 and 3.26, the symbol [ ] indicates evaluation 

at time t^ = t - T, 


[ r * [ ] 


3.27 


t. = t. T 
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3.28 


3.29 


t-.i-T 


3.30 


with T given by Eq. 2.6. In particular 

s T = - T ) 

Note the difference between | V f S T | and j V, S I 

| V, 5 T ) = lyS(*,,i.,2„i-T)l 

/as £S H)\ f3S_ 96 zr\ z ,('9$ LL \ z 

Ux, Utf, at, 3t,9r,/ 

|ys| T = I vis/x., 4l> a., t.;| t t T 

©••(fHSlt,,, 

Finally, indicates the surface defined by the equation 

S T = 3 t -J-J. o 3.31 

Note that 2 T is a surface of the three-dimensional space 
(x r y.^ , z ), which depends parametrically upon x, y, z and 
t. 

Note that, for steady state, the surface 2 coincides 
with the surface 2 . For quasi-steady state, the surface 

^ differs very little from the surface .2 . Thus the 

surface 2 will be called "deformed surface of the body 
and wake". 
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■«rA 


By using Eq. 3.25 and 3.26, Eq. 3.24 reduces to 


4 n £■ (f 



J dS d£ f T J_ I dI T 

^ dk t J 


& fv$ . V fJ.) - .1 — — /'j_)"J T ci) T I 

t L ' 'U' aj- dt, dt, [r p JJ T j 7 pi 


3.32 


- Ill [EF] r S 4 V, 

■a 


V 


Equation 3.32 is the desired generalized subsonic Huygens' 
principle. The reduction of Eq. 3.32 to the classical Huygens’ 
principle and to other well known formulas is considered in 
Appendix B. The more general case in which initial conditions 
are also considered is analyzed in Appendix F. 

Finally, in Appendix C, it is shown that Eq. 3.32 is 
still valid when the control point is on the surface Z if 
the convention is made that E = 1/2 on £ , that is, if Eq. 

3.1 is replaced by: 


E = 0 
= 1/2 
= 1 


inside 
on X 
outside 


1 

X 


3.33 


3.4. Steady Subsonic Flow 

By assuming the time derivatives to be equal to zero, 
Eq. 3.32 reduces to the steady state case 
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4-nEcf = _ & [ys.v,<p - aa 2 - M. H. ] J !_ dt 

z 1 7 9x,J r p RS j 

+ *(-?,)-"*§* ±n )| v J_ 4l 

z r 0X, 0x,' > 'J ' |7S| 


3.34 


- Ill r i «« 

\/ ^ 


Note that Eq. 3.34 differs from Eq. 2.6.10 of Ref. 13 in that, 
in Ref. 13: 


- the convention on the normal is opposite 

- the function represents the volume density of 
a fictitious creation of sources 

- the terms which contain have been neglected* 

(see Ref . 13 , p . 33 ) . 

The correctness of Eq. 3.34 results from the fact that, by 
using the well known Prandtl-Glauert transformation 


X 0 = 

Vl-M* 


VI *•'* 


Eq. 3.34 reduces (in agreement with Eq. B.8) to 

4nEf = - § ft ± ^ 

i t J 0 n. , r . 



v 

where n is the normal to Z. , and 


3.35 


3.36 


*The analysis of the order of magnitude of these terms is 
given in the next section. 
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+ l^f, o' t fo.) + ?0 ') 


3.37 


and, having used the relation 



l 

r 0 I V,SI 


3.38 


3.5. Generalized Supersonic Huygens' Principle 

In this subsection, Eq. 3.19 is specialized to the super- 
sonic case for which the Green function is given by Eq. 2.49 
which may be rewritten as 


with 



3.39 


3.40 


where T + and T_ are given by Eq. 2.48. 

The analysis for the supersonic case is obtained from 
the subsonic one by replacing each term with two, the first 
evaluated at time t^ = t - T + and the second evaluated at 
time t^ = t - T_ . The final result is obtained by modifying 
Eq. 3.32 to yield 


77 £ (P = - # 

1 z % 


? - _L 

d$ 

il 1 

dt,\ 

I T+ ^ 

1 dl T * 

a L 

ad 

dt, 

r 

B 

!v,S T *| 

& 1 

\s. 7, if 

) _ 

M. 

df 

T - , 

, > 

" JS 

r- T. 

€ 

cU, 

cllt, . 

r B 

|v,S T - 


(Eq. cont'd.) 
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* i 4 4UIJN'- 

-1 $ [v,$.V,T-i 4JL (\ + U UL)] T * (fi T * I 1 dZ T ' 

zt f T# L ^ 1 “ ®*.'J V T z |v ( S T *| 

- £ <ft lv,s . - J_ A + u, — )1 f T " 1 — 1 — ^I T 

- |l([£Fl\ [ef] t -) ± mv, 


3.41 


a«J 

■'s* ”T+ " r_ 

and t- are, in general, 


It may be worth noting that 
two different surfaces. 

For steady state S is independent of t and Eq. 3.41 
simplifies into 

_L dz 


1 n Ed> = _ $ M 2, 2S i£ ] ± 

7 *• 2x > -* r R 


r B tfs| 


+ §[*■*(#-«' ff, k ( i)j -pi* 


3.42 


Note that Eq. 3.42 differs formally from Eq. 3.34 only for 
the factor 2n instead of 4n on the left hand side. This 
is a consequence of the well known fact that the supersonic 
doublet is equal to twice the subsonic one: this is due to 

the fact that, as mentioned in Subsection 2.4, the supersonic 
inverse transformation is a one-to-two transformation. 

It should be noted that, according to Eq. 2.40, a point 
can have influence only on the points located inside the 
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"downstream Mach cone" defined by 


* - x, > B 

Thus, it should be understood that the integration is extended 
to the part of the surface for which Eq. 3.43 is satisfied. 

More precisely, G must be considered as a generalized 
function (or distribution according to Ref. 7) defined as 

<3 = - (S^-S T ) *-«, (*,-#]* 

= 0 «-«. 3-44 

This implies that the supersonic doublet cannot be integrated 
in the ordinary way, but must be integrated as the distribu- 
tion theory indicates, that is, the Hadamard finite part of 
7 

the integral must be considered (see Appendix H) . 
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SECTION 4 


SIMPLIFIED FORMULATIONS 


4.1. Introduction 

The formulae derived in Section 3, for subsonic flow, 

Eq. 3.32, and supersonic flow, Eq. 3.48, are very general, 
in that only the hypothesis of potential flow was assumed. 
Usually in practical applications, it is possible (and con- 
venient) to introduce more restrictive assumptions, which 
allow considerable simplification of the above mentioned 
formulae . 

In particular, one can make the assumption of small- 
perturbation flow, which is discussed in Subsection 4.2. 

Furthermore, in most of the practical applications, the 
surface is almost fixed in space. Thus the case of almost 
fixed surface in small perturbation flow is discussed in 
Subsection 4.3. 

Finally, in flutter problems, the surface of the body 
is vibrating with exponentially damped or growing oscillatory 
motion for which the boundary condition can be written as 

H ■ «. ■ 3 - 

where 


s - + tcj 4.2 

is called the complex frequency. This problem is considered 
in Subsection 4.4. Note that, in particular for M = 0, one 
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obtains the harmonic motion. 


It is important to note that the results obtained in 

this section do not require the general formulation of 

Section 3. In other words, using the procedure employed 

in Section 3 under the restrictive hypotheses considered 

here, yields the same results through a much simpler 

derivation. However, it is felt that going through the 

complicated general formulation and then simplifying under 

restrictive hypotheses gives a better understanding of 

the error introduced with the simplification. 

4.2 Subsonic and Supersonic Small Perturbation Flow 

As mentioned in Subsection 4.1, the results obtained 

thus far can be simplified considerably if the hypothesis 

of small perturbation is made. In the following, it is 

assumed that M < 1 (subsonic) or M > 1 (supersonic) but 

2 

not M » 1 (hypersonic) nor |M - 1| « 1 (transonic). 

Consider the boundary condition, which can be rewritten as 


j_ dS i 4>3 

tU |FSJ 

Assuming that (where L is a characteristic length of the 

problem) is small, say of order £ « 1, 

l -- 0(e) 4.4 

L 


is equivalent to assume that the right hand side of Eq. 
4.3 is small, also of order £ 


dS 


V. I VS/ 


oil I 


4.5 
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Thus Eq . 4.3 may be rewritten as 



4.6 


where Q n ~ 0( U . Hence, it is convenient to find an 
asymptotic solution for cf> , as 

<f = 0 (Z**' ) 4 . 7 


In particular, by limiting the analysis to N = 1, one may 
write 

Ip = £ 3 + OfC 1 ) 

I 1 4.8 


It may be noted that all the nonlinear terms are of order 
of £ 2 or higher and thus one can write 

F » F + 0(Z 5 ) 4.9 


By combining Eqs . 4.6, 4.8 and 4.9 with Eq. 3.32 yields 



dL T 


+ 



4.10 





dz 1 + o (rj 


Neglecting terms of order £ z in Eq. 4.10 and returning to the 
original variable yields 
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Equation 4.11 represents the small-perturbation subsonic 
Huygens' principle. 

Similarly, by combining Eqs. 4.6, 4.8 and 4.9 with 
Eq. 3.41, one obtains, for supersonic flow 

4 -fl £ <p(x r *f, Z, t) 



r 8 J lyS 7 'l 


l i t - i7,sr- 

8 J lMS T -j 


*fz. 


T. 





T+ 


r 

B 


iv.sT- jz'- 

K sT ’l 



IV,S[ T - ^Z 7 - 

|V,S'- 
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It should be noted that Eqs. 4.11 and 4.12 are not the exact 


solutions of the linear problem but an approximate (small 

perturbation) solution of the nonlinear problem. For, not 

only the nonlinear terms have been dropped, but the linear 

terms containing dS have been eliminated as well. It may 

cCt 

be noted that one of the reasons to carry the effect of the 
nonlinear terms in the general analysis is to show that the 
order of magnitude of the terms which contain is °f the 

same order of magnitude of the nonlinear terms. Thus, these 
terms can be consistently eliminated if the nonlinear effects 
are neglected. 

It should be noted however, that once the value of cj> 
has been obtained by using Eqs. 4.11 and 4.12, then the effect 
of the neglected terms can be obtained by studying the equa- 
tion of terms of order t. 2 ’ in Eq. 4.10, which has ' cf> z 
as unknown (second order term of Eq. 4.7). It may be worth 
mentioning that this equation contains only known nonlinear 
terms and thus is linear with respect to the unknown 
It may be expected however, that the solution obtained in 
this way is not uniformly valid. Singular perturbation methods, 
which yield uniformly valid solutions, are now under consider- 
ation. 
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4.3. Quasi-fixed Surface in Small-Perturbation Flow 


In many practical applications, the surface of the body 
is almost fixed with respect to the frame of reference. This 
is the case, for instance, of small elastic vibration of the 
wings of an airplane.* In this case, the surface may be 
considered to be fixed in space, although the time derivative 
cannot be neglected. Mathematically speaking, it is assumed 
that 

& = S ^ S(a, i) 4.13 


and 


\7S| r - |VS r | = I vsj 


although 


VS 

0t 




in the boundary conditions. 

Under these hypotheses, Eq. 4.11 reduces to 

-|T 


4nE<f - - J> [§£] rp 

- f Wlty* 

- * m 


’ dZ 


T ± II dz 

> 3*1, 


4.14 


4.15 


4.16 


♦Important exceptions are, for instance, the problems of 
helicopter blades and SDinning missiles. 
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where is the surface defined by the equation S = 0 and 

is given by Eq . 1.13. Similarly, for supersonic case, 
Eq. 4.12 yields 

T. 


4n £f(x, 7, 2, t; = - 


-f 



_L dl 


f(W U <l€)£;(i) dl 

M l UL JZ 


4.17 


9£ 

St 


dt J 


r B So, 


with 9<|) /0d, given by Eq. 1.13. 

4.4. Complex-Exponential Flow 

As mentioned in Subsection 4.1, in flutter problems, 
the motion of the surface is exponentially damped or growing 
vibration. In this case, the surface of the body can be 
written as 

, ( A , , St 

2 = ^ (MJ ♦ L *, (Mj e 4.18 

or, in general 

where S s represents the steady state geometry, whereas 

A 

gives the unsteady contribution. For <f « 1, the 
hypotheses assumed to derive Eqs. 4.16 and 4.17 are valid. 
Hence, Eqs. 4.16 and 4.17 can be used with the surface 2 
described by S g = 0 and boundary conditions given by Eqs. 


4.19 


1.13 and 1.14, with |^5 ( replaced by VS 


, that is. 
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9cf 

2>r) 


where 





4.20 



4.21 


4.22 


with 


"■'J 




4.23 


Since Eqs. 4.16 and 4.17 are linear, 
problems can be studied by setting 


with 





the steady and unsteady 

4.24 


4.25 


and separating the two problems, solving them independently 
and finally, superimposing the results. 

Hence, only the unsteady state component is considered 
in the following. For the subsonic case, by combining Eq. 
4.16, 4.22 and 4.25, one obtains 

4 n£<f e~ S t dL 



4.26 
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or 


4.27 


d„Ef = - f 5 S'jl . f <j JL (■*£- 

Z. "a - 

It may be of interest to note the similarity between 
Eq. 4.27 and the solution of Laplace's Equation given by 
Eq . 3.41. 


Similarly, for the supersonic case, Eq. 4.17 yields 

477 t 2> ^ UL 

21 r - 



where, according to Eq. 2.10, 


-sT. _sT_ 

e e 


- s M( x-*i) 
a B 1 

e - 



4.29 


or, for s = icU (harmonic motion). 


. sT + -ST. 
C + «■ 


, 2 e 


- ^ C x - x; j 
* £> l V ; 


Co S 


cor p 

« B 2 


) 


4.30 


Finally, it should be remarked that Eqs. 4.27 and 4.28 
are suitable for comparing this formulation to the lifting 
surface theory (see Section 5) . However, from the numerical 
point of view, it is convenient to use a slightly different 
procedure described in Subsection 6.6. It should be noted 
that the two formulations differ only for terms of the same 
order of the terms neglected in Subsection 4.2. 
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SECTION 5 


LIFTING SURFACE THEORIES 


5.1. Introduction 

As mentioned in Section 1, the methods currently available 

for evaluating the pressure acting on a wing are based on the 

assumption of zero-thickness-wings (lifting surface theories) . 

The lifting surface theories generally used for aeroelastic 
14 

application are those given in Refs. 8, 15 and 16. In this 
section, it is shown how lifting surface theories are related 
to the formulation presented here. It should be noted that, 
for the case considered here, all the hypotheses of Subsection 
4.4 are satisfied. Thus the results obtained there will be 
used in this section. The subsonic flow is considered in Sub- 
section 5.2 and the supersonic in Subsection 5.3. 

5.2. Oscillating Wing in Subsonic Flow 

Consider an oscillating thin wing in subsonic flow. If 
the thickness approaches zero, then Eq. 4.19 with s = iw re- 
duces to 



where A u 



5 


^ f -- 


5.2 
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and 2. 
the wake . 
plane z, 


is the upper side of the surface of the wing and 
By assuming that this surface is contained in the 
= 0, Eq. 5.2 yields 


rsj 


4ncf = 



5.3 


since 

-“° T — , -i'u)7 

2_ _e _ _ fS ) 

~~ Vz \ y } 5 - 

Note that = 0 outside A ^ ; thus it is convenient to 

replace the domain of integration, 21 , with the whole 

plane (in conformity with theory of distribution). Then, 
differentiating with respect to 2 yields 



dS 

This is a relation between the downwash — 3 L - 

92 

Since the downwash is known at Z = 0, Eq. 5.5 as 2 goes to 
zero yields an integral equation relating the downwash with 

rsj 

A tp . However, for practical applications it is convenient 
to have a relation between the downwash and the pressure dis- 
tribution 


and A f 


^ P ~ r 


According to Eqs. 1.17, 4.18 and 5.6, Ap is given by 



5.7 
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with 


* - P.4, 0^<P ’-4 t ) 

Equation 5.8 can be rewritten as 


. , K ^ / ttO* ^ \ 

- - L'i «• ^(e u - ^f) 


coJX ^ 

e u - a <. 


5.8 


5.9 


By integrating Eq. 5.9 with the condition 




O 


one obtains 


-v) 

= 


- x 

u. 




i‘W \ 

u 4 


p- 


u 


A ?(K^)K 


5.10 


5.11 


Combining Eq. 5.5 with 5.11 and integrating by parts, yields 

_ r^/ 


** ' *' id A» 


4-/1 rr = _L [ C e u - 
32 p V* 1 I ) J 

[*t> « 




- ^ - *0 


- i ' u ;*, 

U * 1 

-< U)T 
£ \ 

■ e 5?( 

*, i 

j 


5.12 


P - 1 


a * oO 

\ l| *.(*•"*» Vi )^.^l 




with 


iWI, 


k . e ‘ 


c (O \ 

y - 2 Z tV 




t £_ 

di 2 - 






4 


5.13 


Note that the finite terms of the integration by parts are 
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equal to zero because of Eq. 5.10 and the condition 


k, = O oX x, = * . 5.14 

which is obtained from Eq . 5.13. By using the transformation 
A - X - A, Eq. 5.13 yields a simpler expression for 



x 


5.15 


with 



Finally, as z goes to zero, Eq. 5.12 reduces to 


5.16 



5.17 


with 

k = ) 5.18 

in agreement with results given by Watkins, Rynyan and 

g 

Woolston . An explicit expression for K , given by Eq. 

D.4S , is derived in Appendix D. Note that in Eq. 5.17, 
the integration can be limited to the surface, , of wing, 

since = 0 outside the wing. 

5.3. Oscillating Wing in Supersonic Flow 

For the supersonic case, it is important to note that 
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the integration domain must be restricted to Mach cone defined 
by 


X-*. > s i/Of- if, )* + **• 


5.19 


Thus, the supersonic case can be handled in a way similar to 
the subsonic case if a cutoff or unit function is used as a 
factor : ^ 


H = H (k-a, - b J 


5.20 


where H (0) 


is the Heavyside step function 


H(0) - / 6 1 > o 

= o e 


5.21 


By performing the same type of operations described in 


Subsection 5.2, one obtains 



- T<1 


in agreement with the results given by Watkins and Berman. ^ 

However, if the wings have only supersonic edges, then the 
two sides of the wing become independent and by assuming a 
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"symmetric flow". 



2 * 


one obtains 





in agreement with Refs. 10 and 16. 


5.24 
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SECTION 6 


NUMERICAL FORMULATION 


6.1. Introduction 

In the preceding sections, a general theory of unsteady 
compressible potential aerodynamics is presented. In Section 
1, the problem is formulated. Section 2 deals with the Green 
functions (for subsonic and supersonic linearized equations) 
which are used in Section 3 to derive an equation which re- 
lates the value of the potential at any point in the field, 

the values of d) and on the surface of the body and 

\ 'b* 

the value of on the wake with an additional contri- 

bution of the nonlinear terms. In Section 4, the formulation 
is simplified for (1) small perturbation, (2) almost fixed 
surface, (3) oscillatory motion, whereas in Section 5, it is 
shown how the classical lifting surface theories can be de- 
rived from the general formulation. 

It is obvious that the general formulation presented 
here has no closed-form solution except for a few very 
special cases. Hence, in general, the use of high speed 
computers will be required. Thus, the numerical solution of 
the problem as formulated in Subsection 4.4 (small pertur- 
bation flow around an oscillating wing) is discussed here. 

It should be emphasized that this discussion is given 
in order to focus a few difficulties which may arise during 
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the numerical computation. Furthermore, it may be noted 
that the difficulties of the unsteady flow (Eq. 4.19) are 
similar to the one of the steady flow. Hence, for the 
sake of simplicity, the numerical formulation is discussed 
for the steady incompressible flow. However, the general- 
ization to the unsteady compressible flow is indicated. 

It should be mentioned that the numerical problems 
arise especially on the treatment of a thin wing (see Sub- 
section 6.4). Thus, in the discussion, it will be assumed 
that the body under consideration is a thin wing, although 
the formulation is valid for any body with sharp trailing 
edges (see Subsection 6.3). 

6.2. Integral Equation Formulation; Existence and Uniqueness 
of the Solution . 

For steady incompressible flow, Eq. 4.27 reduces to 
the classical equation 



Z Z. 


In order to analyze the question of uniqueness of the solution, 
the surface Z. is replaced by a smooth surface Z. sur- 
rounding (at very small, but finite, distance: the boundary 

layer thickness for instance) the body and the wake (Fig. 3a). 
The wake is truncated at a very large, but finite, distance 
from the wing. If the control point is on the surface 2. 
the function E assumes the value 1/2 and equation 6.1 
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reduces to 


ZTItP = - <$> 21 1 dt + 

1 J r 
z ' 

If the geometry of the wake is known and if — — is re- 

0 A ( 

placed by the value (see Eq. 1.12) 





6.2 


which assumes on the body and the wake, then Eq. 6.2 is an 
integral equation relating the downwash integral to the 
unknown value of cp on the surface. Note that Eq. 6.2 
gives the solution of the exterior Neumann problem and, 
in this case, the solution of the equation exists and is 
unique* for any smooth (Lyapunov) surface (Ref. 17, pp. 


620-621). 

Next, the surface 2. (surrounding the body and the 
wake) is replaced by the surface 2. , composed of two 

branches, the surface 2^ of the body and the surface 
2.^ of the wake (Fig. 3b). Thus, Eq. 6.2, combined 


♦Note that, for the interior Neumann problem, the solution 
of the equation is not unique, for any arbitrary constant 
can be added to the solution. Physically speaking, one 
might say that, for the exterior problem, this arbitrariness 
is eliminated by the condition (0 = 0 at infinity. 
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with Eq. 6.3, reduces to 


2n 



2-3 


ds i dL 
0x, r |^s/ 



6.4 


ffo- !«)&(?-) ^ 


2 tul 

w 

In Eq. 6.4, 2 ^ is the upper side of the surface of the 

— ♦ 

wake and hence, the normal fl, is understood to be the 
uooer normal. Note that, according to Eq. 1.17, for steady 
flow. 


£ 

0X, 




on the wake 


6.5 


for no .pressure jump is possible through the wake. 

From physical considerations, the solution of Eq. 6.4 
is "very close" to the one of Eq. 6.2. Thus, it will be 
assumed that, if the geometry of the wake is known, the 
solution of Eq. 6.4 exists and is unique. It should be 
emphasized however, that this conclusion is based upon 
physical reasoning. However, this reasoning is question- 
able, as shown by the remarks given in Subsection 6.5. 
Hence, a rigorous mathematical proof of the existence and 
uniqueness of the solution of Eq. 6.4 would be highly 
desirable . 

However, there are still two important questions to 
be considered: first, the geometry of the wake and se- 

cond, the special behaviour of Eq. 6.4 when the thickness 
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of the wing goes to zero. These two questions are discussed 
in Subsectiop 6.3 and 6.5, respectively. 

6.3. The Wake 


As mentioned in the preceding subsection, the surface 

of the wake in Eq. 6.4 is not known. Thus, Eq. 6.4, which 

is satisfied on the body, must be completed by the equation 

on the wake, which says that the velocity on the wake is 

tangent to the surface of the wake. Thus, one obtains two 

coupled integral equations, one on the body and one on the 

wake, with cP unknown on the body and unknown on 

I 3 o 

the wake, whereas ^ is known on the body and - (^ 

is constant along the x-direction on the wake. According 
to the Kutta condition, this constant value is equal to 
the value of at the trailing edge. Given the velo- 

city on the wake, the geometry of the wake is obtained by 
the condition that the velocity is tangent to the wake. 

This approach has been successfully used in Ref. 6 to study 
the transient incompressible flow around a wing after a 
sudden start. However, from a practical point of view, 
this approach is too lengthy and a simplified treatment of 
the contribution of the wake is presented in the following. 

Note first that 
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where 


ell 



r) 


6.7 


is the projected area (into the plane normal to the direction 
r ) and 


dll , 


6.8 


is the solid angle (see Fig. 4). 

Next, consider the wake integral as a sum of M strips 
in the x direction (see Fig. 5). Applying the mean value 
theorem, one obtains (note that A is only a function of 

y) 

I- * II - - Jl 

' <-(u) 

** W 2 W 

=. _ X IJ (<-{„) fl m 


where are the mean values of Al$> for each strip 

L „ , and are the solid angles of the strip X w 

Equation 6.9 shows that any changes of the wake such that 
solid angles _fl are not altered, do not have any influence 

•'I 

on the value of the wake integral I w . 

This suggests that a "reasonable" geometry for the wake can 
be assumed, provided that the values of the associated solid 
angles are not excessively different from the true ones. Hence, it 
is possible (and convenient) to approximate the wake by straight 
vortex-lines, parallel to the direction of the flow, emanating 
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from the trailing edge of the wing. For, geometrical 
considerations show that the solid angles, _/2 m , are 

changed only slightly. With this assumption, the wake 
integral simplifies considerably and its contribution 
reduces to a line integral. For if the trailing edge 
is given by 

* ~ *t * ( y ) 

6.10 

2 - *r 6 (H) 


then the equation of the surface of the wake is given by 

& * Z T6 (H.) -- O 6.11 


and 


k/Z -» 


1^ * -d^(^S-r)_L dx.dij, 


- h/Z 
b/i 


J Af J w \ 


-i/t 


6.12 


where b is the span of the wing and 

<*5 

J = }(? - 2 ) - ^2 r .- (U .U ll f JlL d.K 


w 


H - y \ 




(fr'lf+te,-*)* 


6.13 


/ - / [*» - * 

'\ r /»,--x Te -(‘ u 


with ? ^ = ? te ('-| | ) . In particular, if trailing edge is 

in the plane 2^ = 0 (i.e. 2 T£ ( ij ) = 0) , Eg. 6.13 reduces to 


J 


w 


at 


r . 

*re ~ * 



i 

(>r E - X )^ + (VH) t r 


6.14 
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In conclusion, under the reasonable assumption of cylindrical 
wake (straight vortex-lines) the effect of the wake sim- 
plifies considerably and Eq. 6.4 reduces to 


2 n 


f * 


& 2S i 

dl. 

+ 

J 9k, r 

|VSj 


J./2- 



/ 

J 

vJ 






6.15 


with J w given by Eq . 6.13. 

Finally, an important remark about bodies without 
sharp trailing edge must be made. In the discussion pre- 
sented in this subsection, it was assumed that the wing 
had a sharp trailing edge. Note that the results can be 
easily generalized to the case of general bodies with 
sharp trailing edge. However, for bodies without sharp 
trailing edge, the inviscid flow theory is incapable, in 
general, of predicting the location of the stagnation 
point from which the wake emanates. This can be easily 
seen in the case of rotating cylinder of finite length. 

From the experiments the location of the stagnation point 
depends upon the angular velocity, to . On the other hand, 
the equation of the geometry of the cylinder does not de- 
pend upon oo and thus oo does not even appear in the 
equation of the inviscid flow. 

In the following, it is assumed that the body under 
consideration has a sharp trailing edge. However, an 
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unsteady viscous theory which predicts the location of 
the stagnation point from which the wake emanates is neces- 
sary in order to extend this method to bodies without 
sharp trailing edges. Similar consideration holds for the 
case of detached flow (when this can be approximated by 
a wake emanating from a point different from the sharp 
trailing edge). 

6.4. Numerical Solution of the Integral Equation 

In order to solve Eq. 6.15, various approximate 
techniques are available. For lifting surface theories, 
one of the most successful ones is the collocation method. 

In this method, the unknown function is approximated by 
a linear combination of N prescribed functions with un- 
known coefficients. The functions are generally the first 
N ones of a complete set of functions, each of which satisfy 
the boundary conditions of the problem. The N unknown 
coefficients are determined by solving a linear system of 
N equations obtained by satisfying the integral equation 
at N points (collocation points) . This method is being 
explored for the solution of Eq. 6.15. However, it should 
be noted that, for complex geometries, this approach is 
not feasible, for each geometry requires a different set 
of functions and the more complicated is the geometry, the 
more difficult it is to guess the appropriate set of 
functions . 

Hence, a more flexible approach is desirable. A 
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representation of the unknown function similar to the one 
used in finite elements for structural problems seems to 
be more convenient. In finite elements, the domain of the 
equation is broken into small elements. The function in- 
side the element is expressed in terms of its unknown values 
(and the values of its derivation, eventually) at nodes of 
the element. This general tvoe of representation is now 
under examination. It should be noted that, in finite 
elements, the final eouations are derived from variational 
principles. Here, the same type of representation is used, 
but the algebraic equations are obtained by satisfying 
Eq. 6.15 at prescribed points (control points). 

The most elementary form of this approach is very 
close to the box method and is described in the following. 
Consider the surface Z! divided into small elements 
(.see Fig. 6) 


Z. 




6.16 


By the mean value theorem 


HAM* 


where cp. is an appropriate value within the element 
This suggests that ^ may be approximated by the 
value, (P. , at the center of the box. This yields 
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<&>- # a zb 4,^n^) di - 


z 


6.18 


4. Z (P VJ 
H, H 


t l 


where, in the last term, the index l covers only the boxes 
in contact with the trailing edge and 


II 

1+ 


6.19 



where the upper (lower) 

sign is used for the upper (lower) 

side of the wing. 



By satisfying this 

equation at the centers 

P (k) of the 

boxes, one obtains 



H 



(P = t + 21 

c cp + Z. vr 'if. 
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where 



1 

ii 


* 6.23 


is the distance of the dummy point of integration from 
(k) 

the center P of the element k. Finally, 


u = W 


Ki 


ip-- p 
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WO 
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for the boxes in contact with the trailing edge and 


= o 

otherwise. Equation 6.20 can be rewritten as 


ri 



with 


6.25 


6.26 





6.27 


where S fc . is the Kronecker delta. Eq. 6.26 is the 

equation which yields the solution of the problem. 

6.5. Limiting Behavior for Zero Thickness 

As mentioned above, the formulation described thus 
far becomes singular in the case of zero thickness. This 
is shown clearly by the fact that, for lifting surface 
(in the plane z ±= 0). In this case, Eq. 6.1 reduces to 


In 


t 


!<*.- v s (?) 

-u, 2 ' i0 


6.28 


and 


rH a -=° 


' — li 


6.29 


"x ("'I 

where 2. is the portion of the plane = 0 (upper side) 
which contains the wing and the wake. By adding and sub- 
tracting Eqs. 6.28 and 6.29, one obtains 


1 ♦ i - 


, 0 


6.30 
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and 



6.31 


This implies that (since, as well known, there exists a 
nontrivial solution ) the operator shown in Eq.6.31 

is singular. 

Hence, one can expect that the numerical procedure 
also has a singular behavior. In order to show that this 
is indeed the case, consider a symmetric wing with angle 
of attack and thickness ratio t , and let t go to 
zero. In this case, Eq. 6.21 shows that 


V\ya t> k = 0 6 .32 

T-*o 

In order to simplify the discussion, the numbering of the 
boxes is assumed to be such that the odd (even) numbers 
correspond to boxes in the upper (lower) surface and that 
the box in opposite position to the upper box i, has the 
number i + 1 (see Fig. 7). For simplicity, upper box i and 
lower box i + 1 will be called "opposite boxes". 

With this numbering, it is easy to show that, according 
to Eq . 6.22, 
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In other words, all the coefficients are equal to zero 

except for the ones relating opposite boxes, which assume 
the value -1. Furthermore, the coefficients are equal 

to zero. Hence, Eq . 6.26 in the limit, as x goes to zero, 
reduces to 


i i : o o ; .io o 

f i i 

i l i O O < j o o 

b O > ( I I j o o 

O O • I 1 \ o o 



= 0 
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0 o ■' o o I 

O O | o O ; 


This equation can have nontrivial solution since the deter- 
minant is equal to zero. 

Note that this result implies that zero thickness 
wings (lifting surface theory) are more difficult to deal 
with than finite thickness wings. 

However, this shows also that, by using the method 
proposed here, one may encounter numerical complication 
due to the fact that, for very thin wings, the determinant 
is close to zero and hence, one may encounter strong eli- 
mination of significant figures. This implies that one 
has to be very accurate in the evaluation of the coefficients 
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and . 


It may be of interest to analyze the order of magni- 
tude of the different terms of Eq. 6.26: by writing Eq. 

6.26 as 

Hh K 


t» ? 


6.35 


and noting that 

[(fj * °('j 

l = 0(tJ 

l * ; 


6.36 

6.37 


one obtains 




6.38 


In order to establish the practical limits of the 
applicability of the method, Eq. 6.26 has been solved 
numerically for very small values of T . The results 
are presented in Section 7. 

6.6. Generalization to Unsteady Subsonic Flow 

In this subsection, the formulation presented above 
is generalized to cover steady and unsteady subsonic flow. 
As shown in Appendix C, great simplification is obtained 
if the generalized Prandtl-Glauert transformation 

= t 

Ho 8 W 

*. = z 6 * 39 

to -- 
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is introduced. Following the same procedure used in Sub- 
section C,3.2, by using Eg. 6.39, neglecting nonlinear 
terms. Eg. 3.32 reduces to 


4n£ f = 




_M 11. _ M ds gf 

3 5x ol 2X ol 


T * , dz* 


r ' l v J r 1 


+ )-n$i -Z-fi) 

#*■ 


(f z dll' 


I^S 1 


at, ^ 
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-iT. 




dZ. r * 


r » (Y.S* 


6.40 


where 




6.41 


and 


I M L,. 


6.42 


with 


T. = = M(^-\) + r. 


6.43 
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Next, it is assumed that the surface is almost fixed and 
Eqs. 4.13 to 4.15 can be used, to yield 


4n£<p z _ & f*L _ M asyfrs 

" itj,s| d t,, 


1 dl 0 

r a 


+ 




6.44 


- & 211 1 dl* 

£ 3o Ci fit. 


Finally, neglecting terms of the same order of the nonlinear 

terms (which implies that ®}2. can be replaced by ^ = 0„ ) 

0a, 9n ‘ 

Eq. 6.44 reduces to 

4„E f * -§ Qj- X «£ 




6.45 




0O, r, 0to 


This is the desired equation. It may be noted that M 
appears only in T Q , whereas (2> does not appear at all. 
For steady state, Eq. 6.45 reduces to 


4n£ f - - § + #£, (-r) f ^ 


6.46 
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This shows that the same method used for incompressible 
flow can be used for steady compressible subsonic flow. 
This equation is essentially the same as Eq. 2.6.10 of 
Ref. 13, already discussed in Subsection 3.4. 

Furthermore, for unsteady oscillating flow (as des- 
cribed in Subsection 4.4), combining Eqs. 4.1, 4.18 and 
6.43 and 6.45 yields 


4 71 £ <f - - jj> 
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n r. 
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where 


s. = — 
I s a - 


6.48 


This equation can. be simplified further by introducing the 
functions 


G. - 0. e 

X* 

t - f C 


6.49 


6.50 


-s.r. 


to yield 

A if A 

= - C _e 


I i &(4' 


6.51 


*12. 


This equation is equivalent to Eq. 4.27: the difference 

between the two equations is of the same order of the terms 
which have been neglected. 

Furthermore, Eq. 6.51 shows that., by using the general- 
ized Prandtl-Glauert transformation (Eq. 6.39) the equation 
relating ^ to is completely independent of Mach 

number. Note however, that the contribution of the wake 
depends explicitly upon the Mach number, M: for, using 

Eq. 6.34, Eq. 5.9 reduces to 


p. u : 
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where 


.A 

Ap - Af> e 


i.nx, 


6.53 


which implies that, on the wake, where Ap = 0, 

>o X. S. *oTE 


A 


Ad I e 


ir 


(“l 


6.54 


Note also that the only terms neglected are those in 
the integral which contains the boundary conditions. This 
is important, because, as shown in Subsection 6.5, the 
operator becomes singular when the thickness goes to zero. 
Hence, even small terms may become important when the 
thickness becomes small. 

In conclusion, Eq. 6.5 is more suitable than Eq. 4.27 
from the numerical point of view. The procedure used to 
solve Eq. 6.1 can be used with minor obvious modification 
to solve Eq. 6.51. The only complication arises from the 
contribution of the wake, which can be treated as described 
in Appendix D. 


65 



SECTION 7 


NUMERICAL RESULTS 


7.1. Introduction 

As shown in Subsection 6.5, in the case of very thin 
wings, one may expect to encounter strong elimination of 
figures. In order to establish the limits of applicability 
of the proposed method, Eq . 6.26 has been solved numerically 
for a simple case for which numerical results are available: 
rectangular wing in steady subsonic flow. For the sake of 
simplicity, the procedure is described for incompressible 
flow only, since, by using the Prandtl-Glauert transformation, 
the steady compressible flow reduces to the incompressible 
one (Eq. 6.46 in Subsection 6.6). 

7.2. The Geometry of the Wing 

Consider a rectangular symmetric wing with thickness, 
h, given by 

K = x c 3|I fg (l-| ) fj -Y 7.1 


with 


€ - * /c 


7.2 


ij = 2 Cj / b 

where c is the chord and b is the span, x and y are the 
cartesian coordinates of the planform at zero angle of attack 
(see Fig. 8) and, finally, 
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I 


7.3 


T - 



U 


T-O 


is the thickness ratio. 

Equations 7.1 and 7.2 can be used to give the geometry 
of the wing in parametric form (parameters § and ) , at 

zero angle of attack, as 
* = c | 

i ■ 1 7 

l = 4 t ci5 (i-|) 

4 

= i h. 

2 

where the upper (lower) sign holds for the upper (lower) 
surface of the wing. 

If the angle of the attack, <*. , is different from 

zero, the geometry of the wing is given by (See Fig. 9) 

X c X CoS bi +2 Silo L 


1- 1 


7.5 


For small values of z and 
approximated as 

x = x 

1 - i 

2 = i - K U 


Eq. 7.5 can be 


7.6 
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Hence, the surface ^ 
which implies 


can be written as ££-h.±xtx.-0 

2 


9-S . [ +• rJ 

0x ~ 2 a* 


dS 


J jH/ 

2 a j 


2s 

2? 


± I 


7.7 


Equations 7.6 and 7.7 fully describe the geometry of the 
wing and enable one to evaluate the coefficients C k . and 

7.3. The Numerical Procedure 

The numerical procedure used to evaluate the coefficients 
of the equation is briefly described in the following. First, 
note that the wing is symmetric with respect to the plane 
y = 0. Hence, Eq. 6.26 can be rewritten as 

hi fj 



A 

where N is the total number of boxes on the right hand part 

of the winq and c is the influence of two boxes (in 

M 

symmetric position with respect to the plane y = 0) on a 
point l on the right hand part of the wing, 
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with 


7.9 
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where 

and all the other quantities are evaluated on the right hand 

A 

part of the wing; similar expressions hold for w . 

pi 

Second, note that 3 S/2x, , and 3S/2tj ( are 

infinite at the leading edge and the tip of the wing, 
respectively. Hence, it is convenient to use a nonuniform 
mesh for the definition of the boxes (smaller boxes in the 
neighborhood of the leading edge and tip, larger boxes in 
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the neighborhood of the trailing edge and root) . This is 
accomplished by introducing the transformation 

§ = X* ( O ± X ^ i ) 


Jj - i-(/-v ) 2 [o ± y ^ i ) 


7.12 


and using boxes of constant sizes A X 
X, Y : 

^ X = I/MX 
AY = 1 / M 


AY in the plane 


7.13 


where NX and NY are the number of boxes in direction X and 
Y respectively. In other words, the center of the box 
(m, n) is given by (see Fig. 10) 


whereas. 


with 
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( *»- i) A * 

M : lj • 

■ , NX 

y to 
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ytrt) _ y(c I _ AV = (,.|jAV 

•» rt 2 


+ 


2V _ n AY 


Note that, for each couple of values m and n, there are two 
boxes, one on the upper and one on the lower side of the wing: 
hence, the total number of boxes on the right hand part of 


the wing is 


KJ = 2 NX NY 


Third, the coefficients bp are evaluated as 


■ 


_ ZL & ±h 

Y~> U Zn \ r r [y: 


7.17 


7.18 


where are the same boxes used for the evaluation of 

Cpg. This procedure is particularly convenient because it 
yields very accurate results and is very little time con- 
suming, since most of the operations required are needed 
anyway for the evaluation of the coefficients 

Fourth, it should be noted that the emphasis here is on 
very thin wings. Hence, it is feasible to approximate each 
surface element with its tangent plane at the center 

of the element; the boundary of the element is still given 
by Eq. 7.15. In this case, the integrals can be evaluated 
analytically. This is shown in Appendix E for the more 
general case of trapezoidal element, which is needed, for 
instance, in the case of swept or delta wings. 
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Finally*, since only very thin wings are considered here, 
the pressure coefficient, given by Eq. 1.18, is evaluated 
as 


(* _ _ z - - Z 2 x 

i 0 K 


9x 


7.19 


- I 
tUx 

where /dx is evaluated by central finite differences. 

Note that [d x is finite whereas 3^/dx is infinite 
at the leading edge. This is one of the advantages of using 
the transformation given by Eq. 7.12. 

7.4. Numerical Results 

In order to compare this method with experimental and 
lifting surface results, the rectangular wing considered in 
Refs. 18 and 19, for which 


* = 5 ' 


7.20 


b /c = 3 
H = . 24 

is investigated here. Analysis of the thickness effect (pre- 
sented in Subsection 7.4.1) shows that the solution obtained 
by employing a thickness ratio z = .001 is a good repre- 
sentation of the zero thickness solution. Furthermore, 
analysis of the convergence (presented in Subsection 7.4.2) 

A 

shows that using NX = NY = 7 (that is N = 98) is sufficient 
for the convergence. The results obtained with Z= .001 
and NX = NY = 7 are shown in Fig. 11 and 12 where the 
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distribution of the potential over the wing and the lift 

distribution 


C - _ Zic =. c , - c* 7.21 

respectively, are presented in threedimensional form. 

It may be noted that the diagram of A ^ is flat in 
the neighborhood of the root of the wing and the trailing 
edge (more precisely QA^/d^ = 0 at the root and 0A*^ /&/■ = 0 

at the trailing edge). Similarly, the diagram of c is flat 
in the neighborhood of the root (more precisely /®H ~ ® 

at the root). Hence, the values of A and c^ at the 
center of the boxes in contact with the root (root boxes 
values) and the value of A ^ at the trailing edge boxes 
will be considered in the following in order to discuss 
the effect of the thickness and the convergence. 

7.4.1. Thickness Effect 

* 

In order to analyze the thickness effect, the problem 

has been solved for four values of the thickness ratio, 

T = .1, .01, .001 and .0001 respectively. In all these 

cases, the number of boxes in both x and y direction is 

NX = NY = 4. Hence, the total number of boxes (for upper 

- . A 

and lower side of the right half of the wing) is N = 32 
(i.e., Eq. 6.26 is a system of 32 equations and 32 un- 
knowns). For the value T = .001, no message indicating 
strong elimination of figures was given, whereas, for the 

*This is not an analysis of the thickness effect on the real 
solution (since that would depend upon nonlinear terms) but only 
the examination of the effect of decreasing thickness on the 
numerical process. 
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value I 


= .0001, a message indicating an elimination of 
significant figures higher than the prescribed tolerance at 
the 19th step was obtained.* 

Hence, only the cases T = .1, .01 and .001 are pre- 
sented here. The values of at the centers of the 

trailing edges boxes and the root boxes are shown in Figs. 
13 and 14, respectively, whereas, the lift distribution 
C, = - Ac at the root boxes, is shown in Fig.- 15. The 
results indicate that the solution converges to a zero- 
thickness solution and that the solution for t= .001 is 
a good approximation for the zero thickness solution. 

Note that, according to Eq. 6.6 





since the point from which the solid angle is evaluated 
is on the surface 21 . This equation is poorly satisfied 

on the leading edge and the tip where the approximation of 
the surface element with its tangent plane is poorer. The 


*For the solution of Eq. 6.26, the standard IBM SUBROUTINE 
GELG has been used. The value of the tolerance (which is 
compared to the ratio between the pivot at the n-th step and 
the initial step) was chosen to be TOL = .001. 
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poorest values of Zj C at "tip or leading edge boxes", are 

^ PP 

given in Table 7.1, column 1, whereas the poorer value for 
the "internal boxes" (not at the leading edge nor at the tip) 
are given in column 2 . 


TABLE 7 . 1 


i r 

l 

2 

.i 

.76946 

.98986 

.01 

.96247 

.99899 

.001 

.99625 

.99990 

. 0001 

.99960 

.99999 


Table 7.1 indicates that the approximation of the surface 
elements with its tangent plane is not satisfactory for the 
"tip or leading edge boxes" for the case T = .1. A more 
sophisticated analysis, which evaluates the error (difference 
between integral on the tangent plane and integral on the 
real surface element) by Gaussian numerical quadrature 
formulae is now being analyzed.* 


*Note that the use of X and Y See Eq . 7.12) as variables of 
integration eliminates the singularity of the integrands at 
the leading edge and the tip. 
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7.4.2. Convergence 

In order to study the convergence of the solutions, the 
case T = .001 was solved for NX = NY = 4, 5, 6 and 7, res- 
pectively. The value of at the root and trailing edge 

boxes are shown in Fig. 16 and 17 respectively, whereas, the 
values of the lift distribution c , = - Ac at the root boxes 
are shown in Fig. 18. 

The results show that the solution is convergent very 
fast and that the case 4 x 4 is sufficient for an accurate 
analysis. The computer time employed on the IBM 360/50 
available at the Boston University Computing Center, are 
given in Table 7.2. 


TABLE 7.2 


Number of 

Computing Time 

Boxes 

Sec . 

4 x 

4x2 

22.2 

5 x 

5x2 

60.4 

6 x 

6x2 

129.9 

7 x 

7x2 

259.9 

Comparison with Existing Results 


In order to evaluate the accuracy of the method, the 
results shown in Fig. 12 are compared to the one obtained 
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in Ref. 19. 


For convenience, three vertical sections of the three- 
dimensional diagram presented in Fig. 12, are plotted in 
Fig. 19. The three sections correspond to values of = .5, 
.7 and .9 respectively. 

In order to evaluate the comparison, the following 
factors should be emphasized. First, this test case was 
considered in order to verify the applicability of the 
method in the worst possible conditions (very thin wings 
with thickness ratio 1 = 1/1000). Second, the numerical 
procedure was chosen for its flexibility (i.e., possibility 
of applying it to very general geometries) and not for its 
accuracy. Furthermore, it should be noted that the com- 
parison should not be made with the experiments, but rather 
with the lifting surface theory, since the thickness ratio, 

"C = .001 is considered here, is very small. Finally, it 
may be concluded that the results obtained here are in 
surprisingly excellent agreement with the ones presented 
in Ref. 19. 

Note that the case T = .1 (which represents a realistic 
value of the thickness ratio) is only partially satisfactory, 
because Eq. 7.22 is poorly satisfied. Hence, a comparison 
with the experiments is not attempted here. As mentioned 
above, a more accurate procedure to evaluate the coefficients 
c and b (in this procedure, the surface elements are not 

pq P 


77 



necessarily planar with straight boundaries) is now under 
* 

investigation. Preliminary results are in good agreement 
with the ones presented here (Ref. 20). 

On the opposite end of the range of the thickness 
ratio (the case X = .0001) strong elimination of signi- 
ficant figure was obtained (see Subsection 7.4.1). However, 
note that this value of the thickness ratio is much too 
small to be of any practical interest. Furthermore, despite 
the elimination of figures, the results obtained were very 
close to the ones for r = .001. In conclusion, there is 
no limitation of the method (at least for cases of practi- 
cal interest) due to the singular limiting behavior (for 
zero thickness) described in Subsection 6.5. It should 
be remarked again that finite thickness wings can be treated 
in a simpler fashion than zero-thickness wings. 

Finally, it should be emphasized that, once the values 
are known, the potential c|> and the pressure coefficient 
Cp can be evaluated at any point of the field. For, by 
using the same procedure applied to derive Eq. 6.18 (with 
E = 1 instead of E = 1/2), Eq. 6.1 can be approximated as 

u 
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4T\ 
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dL 


- -1 2 <p ff f fi) dU 
■4-n f.i I' Jj 
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Similarly 


<L - - ZH - - 4 H 2 (±)4Z 

^ ~ zn jf 0n, 0x W 


7.24 
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- - ?-<p [[ 2L JL/jJdZ . Z (P Oak 
Zn '• Jl 2>r\. &x\r/ ^ *=• T« 2< 


0tl, 

2. 

Note that these expressions can be evaluated by using the 
same trapezoidal elements described in Appendix E: the 

coefficients to be used in Eq. 7.24 are simply the de- 
rivative with respect to x of the coefficient given in 
Appendix E. Furthermore, if the point is on the surface 
T , Eq. 7.23 and 7.24 are still valid if coefficients 
are evaluated in the limit sense and the value of E is 
maintained equal to one. 

Finally, note that the lift and the moment coefficients 


can be evaluated as 

W/i «• 


Mz. 


C CTnI ^ J c t •- £ / A fr 6 ^ 

C *kJ f c t '• &/ •H, j A t. 

-t>/2 ; .*/i i 

w/z. ^ 

: 4_] ( [*f *■] - j -aij ) J*,) cbf t 
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-- ±] C ^r S <L i, 

e*i , i 

Similar relations can be used for evaluating the 
generalized forces. 

- 


7.25 


7.26 
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SECTION 8 


DISCUSSION 

8.1. General Comments 

As mentioned in Subsection 1.1, lifting surface 
theories have two main disadvantages. First, they are 
complicated from the numerical point of view (singularity 
of the integrands and, especially for unsteady state, 
complication of the kernel function). Second, they can- 
not be easily generalized to study complex configurations 
and motions. The formulation presented here is an attempt 
to reduce these disadvantages. The numerical simplicity 
of the proposed formulation, in comparison to the steady 
and (especially) unsteady lifting surface theory, is 
apparent from Sections 6 and 7 . Note that a sufficiently 
accurate pressure evaluation requires only 22 seconds of 
computing time on an IBM 360/50 (see Subsection 7.4.2). 

The second point, applicability to complex configurations 
and motions, is discussed in the next subsection. 

8.2. Applicability of the Method 

The main advantage of the method proposed here is 
the fact that it can be used to solve a large variety of 
problems, for which the lifting surface theories can be 
used only in a very unsatisfactory way. It should be 
noted that the method, although classic, (as shown in 
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Section 5, the lifting surface theories can be obtained 
as a limiting case of the present formulation when the 
thickness goes to zero) has never been developed to its 
full generality. 

Note also, that the method is valid for both subsonic 
and supersonic flow. The procedure described for the sub- 
sonic flow can be used for the supersonic one. However, 
for supersonic trailing edges, the procedure is simpler, 
since the wake has no contribution on the body. 

In order to appreciate the extent of the applicability 
of the method, it might be convenient to consider a few 
typical examples. 

A simple, but interesting application of this method 

is the evaluation of the aerodynamic forces acting on wings 

of finite thickness. The importance of thickness effects 

in the prediction of flutter boundaries has been shown by 

21 

Yates and Bland. It should be noted that, in the pro- 

posed formulation, the nonlinear effect of the thickness might 
be included in a very systematic and natural way, by using 
singular perturbation methods mentioned in Section 4. 

More importantly, the method can be used to solve 
more complicated problems, since it is formulated for 
arbitrary geometries and motions. The arbitrariness of 
the geometry implies that even complete configurations 
(wing, body, tail) can be studied, including wing-tail 
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wake interference effects. Even more detailed configurations, 
like tanks at the wing tips can be analyzed without any in- 
crease of difficulty. 

On the other hand, the arbitrariness of the motion 
implies that problems which cannot be examined with lifting- 
surface theory, can now be solved, although the solution is 
rather complicated. Typical examples are: curved tra- 

jectories, accelerated motion, and roll (which is particularly 
important for practical applications) . Other problems 
like gust response (and indicial motions in general) can 
be solved, in a relatively easy way, by making use of the 
generalized formulation, derived in Appendix F , which in- 
cludes the effects of arbitrary initial conditions. It may 
be noted that damped oscillatory motion, which is important 
in predicting the degree of stability of linear systems, 
as well as arbitrary periodic motions do not offer more 
difficulty than the simply harmonic motion. 

A particularly interesting problem is the evaluation 
of the aerodynamic pressure on the blades of a helicopter 
in forward flight in which both geometry and motion are 
extremely complicated. These examples show the appli- 
cability of the proposed formulation in solving problems 
involving complicated configurations having arbitrary 
motions . 

It is of interest to mention that, by making use of 
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the method of images, described in Ref. 11, the present 
formulation can easily be extended to study the motion in 
the presence of an .infinite rigid planar wall which acts 
like a plane of symmetry: this corresponds to the practical 

problem of motion in the vicinity of the ground (ground 
effect). Also, using the method of images, the problem 
of flow inside a rigid circular duct can be studied, as 
well as the slightly more complicated flow inside a duct 
of rectangular section. 11 These cases correspond to the 
important problem of accounting for the effect of the 
walls of a wind tunnel. In summary, with the method of 
images, the proposed method can be extended to study the 
ground effects and the effect of the wind-tunnel walls on 
the experimental results. 

It may be worth noting again that the formulation 
reduces to appropriate formulae in the particular cases 
considered in Appendix A. Also, the lifting surface 
theories can be derived as a limiting case of this for- 
mulation. 
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SECTION 9 


CONCLUDING REMARKS . 

A general theory of steady and unsteady, compressible 
aerodynamic flow around a lifting body having arbitrary shape 
and motion has been developed. The theory is based upon the 
classical Green function method. This yields an integral-dif fere 
tial representation of the velocity potential. For the important 
practical case of small perturbation, if the control point is on 
the surface of the body, the representation reduces to an in- 
tegral differential equation relating the potential on the 
surface to its normal derivative. In particular, for small 
harmonic oscillations around a rest configuration, one obtains 
a two-dimensional Fredholm integral equation of second type. 

This formulation reduces properly to lifting surface theories 
and other classical results. The question of uniqueness was 
examined and it was found that, for thin wings, although the 
operator becomes singular as the thickness approaches zero, good 
numerical results can be obtained even for thickness ratio, f = 
. 001 . 

In conclusion, the formulation developed here can be used 
for efficient numerical solution of a large variety of problems 
for which no satisfactory methods are available. Hence, the 
method should be more convenient than existing ones, even for 
very simple problems, as the rectangular wing in steady subsonic 
flow presented here. 
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Fig. 11 


Potential distribution for a rectangular wing with aspect ratio b/c = 3, 
angle of attack OC = 5° and Mach number M = .24 



Fiq 12 Lift coefficient distribution for a rectangular wing with aspect ratio 
b/c = 3, angle of attack oC = 5° and Mach number M = .24 


Fig. 13 Effect of decreasing thickness on the numerical scheme: Potential 

difference at the trailing edge boxes for NX = NY = 4 and different 
values of the thickness ratio, ~c , for a rectangular wing with 
aspect ratio b/c = 3, angle of attack oC = 5° and Mach number M = . 





Fig. 14 Effect of decreasing thickness on the numerical scheme: Potential 

difference at the root boxes for NX = NY = 4 and different values 
of the thickness ratio, T , for a rectangular wing with aspect ratio 
b/c = 3, angle of attack oc = 5° and Mach number M = .24 
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Fig. 15 Effect of decreasing thickness on the numerical scheme: Lifting 

pressure coefficient at the root boxes for NX = NY = 4 and different 
values of the thickness ratio, t , for a rectangular wing with as- 
pect ratio b/c = 3 , angle of attack oc = 5° and Mach number M = .24 
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Potential difference at trailing edge boxes for thickness ratio 
T = .001 and different values of NX and NY for a rectangular wing 
with aspect ratio b/c = 3, angle of attack = 5° and Mach number 
M = .24 


Fig. 16 



Fig. 17 Potential difference at the root boxes for thickness ratio T = . 

and different values of NX and NY for a rectangular wing with as 
pect ratio b/c = 3, angle of attack OC = 5° and Mach number M = 



104 



Fig. 18 Lifting pressure coefficient at root boxes for thickness ratio 

' r = .001 and different values of NX and NY for a rectangular wing 
with aspect ratio b/c = 3, angle of attack oc = 5° and Mach number 
M = .24 







APPENDIX A 


TWO FUNDAMENTAL FORMULAE 


In Subsection 3.3, the formulae 


[[ I E AS) Ai 
JjjJ 2L r 1 

- Mf 

and 

flf fjM 'zwJt, 

JfJi 


--§> f T jO£T 

s'" ivs t | 



A.l 


A. 2 


were used (see Eqs . 3.25 and 3.26) to derive the generalized 

Huygens' Principle. In this Appendix, it is shown that A.l 

and A. 2 are valid for any "good function", as defined by 
. 22 

Lighthill. The concept used in the following is very 
similar to the one introduced in Ref. 22. 

In order to prove Eq. A.l, it is convenient to assume 
the "surface distribution", Oj , as the limiting case 
of volume distribution in a thin layer of the four- 
dimensional space (with constant infinitesimal thickness, 
t* ) surrounding the surface S(x, y, z, t) =0. This 
is equivalent to the classical procedure of defining a 
surface integral (which is connected to the distribution 
, see Eq. 3.7) as a limiting case of a volume in- 
tegral. Thus define 
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A. 3 


where the "layer function" L (similar to the function E 
defined- by Eq. 3.1 is defined by 

L = 1 inside the layer 

A. 4 

= 0 outside the layer 

More precisely, the function L has the value one for the 


fill 

points (x^, y^, z^, t^) such that 


1 

as | 

[ X ' - X. 1 

i + 3S, 

< u 1 

I + 2S , ) + 2S /. I_- J £* 

|0,s| 

3x, 

V 1 ij 


Li *ij 

1 diy* *9t, t,] < t a. s 


In Eq. A. 5, all the derivatives of S are evaluated on 

_ iiii 

the surface Z. • Furthermore (x^ , y^ , , t^) is a point 

in the direction of the normal, V , to the surface £ at 
the point (x^, y^, z ^ , t^) . Note that, by definition, 

3 (*,, H, •> ) - 0 A -6 


Combining Eqs. A. 3 and A. 5 yields 


- 06 


I 1 JJJJ I L s r W'M, 


A. 7a 


Performing the time integration yields 

4 5 z ^ t m'A -- i Hi r l t d\>[ 




A.7 b 
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The integrand is different from zero only for the points 

of the layer of the three-dimensional space for which 
T 

L = 1. More precisely, the integrand is different from 
zero for the points for which Eq. A. 5 is satisfied with 

I I 

t^ = t - T and t^ = t - T; that is, for the points for 


which 



Note that the thickness of the layer is infinitesimal, 
and thus neglecting higher order terms, T - T can be 
expressed as 


t'-t = 

2 a, a if, v c, 


A. 9 


where the derivatives are evaluated on the surface 
Combining Eqs. A. 8 and A. 9 yields 



ras 3T I 1 /**.* _ 11 ?L lVu/.tj ) 
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T 

where S is given by Eq. 3.28. Thus, the layer in which 
the integrand is different from zero, is a layer surrounding 

Z t t 

of equation S =0; the thickness of the 

layer is given by 


£ 




A. 12 


Hence, performing the integration through the thickness 
and taking the limit, Eq. A. 7 reduces to the desired Eq. 


A. 1 . 

Next, in order to prove Eq. A. 2, consider 
as a limit of the incremental ratio. Then 



- lirn 
At -»o 



By using Eq. A.l and performing the limit operation, Eq. 
A. 13 reduces to the desired Eq. A. 2. 
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APPENDIX B 


REDUCTION TO ELEMENTARY CASES 


In this Appendix, it is shown how classical results 
can be obtained as particular cases of Eq. 3.32. 

B.l. Huygens' Principle 

The problem related to the Huygens' principle is a 
particular case of the one considered here, the differences 
being that 

1) The velocity of the undisturbed flow is equal to 
zero 

2) The surface is assumed fixed in time 

3) The nonlinear terms are zero 
Mathematically speaking, these differences can be ex- 
pressed as 

U. - O (H.o) 

0S/0k= O B * 1 

F =o 

Note that, according to Eq. B.l 


^ (J ££ _ q 

dik " * ** 

Combining Eqs. 3.32, B.l and B.2 yields 


B.2 


Z 


— L_ dl 

iVl 


B. 3 
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- # [f?f 


iv.si 

dl 
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4Z 

lv,s 
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where 




and 


r = 


(<- if ,) 1 + ((/ -^,) 4 + (2 '*/)* 


B. 4 


B . 5 


Note that 

0, . &- 

|y s i 

is the outward normal, directed from the region E = 0 
to the region E = 1. Thus, Eq. B.3 reduces to 




B. 7 


which is the well known Kirchhoff formula which is the 
mathematical expression of Huygens’ principle (see Refs. 

11 and 12 in which the opposite convention on the definition 
of the normal is used) . 
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B.2. Integral Representation of Solution of Poisson's 


Equation 

The differences between Poisson's equation and the 
case considered here are: 


1) No time dependence 

2) M = 0 


12 

Then Eq. 3.32 reduces to the classical result 



V 


In particular, for F = 0 (Laplace's equation) Eq. 
reduces to the well known formula 


B. 8 


B.8 



B. 9 


B.3. Poisson's Formula 

In Appendix F, it is shown how Eq. 3.32 can be used 
in order to derive the contribution of the initial con- 
ditions. Here it is shown that for M = 0, the contribution 
of the initial conditions reduces to the well known Poisson 
formula, that is, the solution of the wave equation with 
given initial conditions. The contribution of the ini- 
tial conditions is derived in Appendix F and is given by 

Eq* F.2. For U= 0, one obtains r, = r, T= — r and 

* Pa 
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Eq. F.2 reduces to 



where IL* = t - T is the spherical surface of equation 

r - B. 11 

By using Eq. B.ll, Eq. B.10 may be rewritten as 



■L 

4/1 




V0t. 4. 


4-n 2t 



B. 12 


where d SL = dl-*/r l . Equation B.12 is the well known 
12 

Poisson formula 
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APPENDIX C 


THE VALUE OF THE FUNCTION E ON THE SURFACE 


C.l. Introduction 

Consider Eq. 3.32: neglecting the nonlinear terms 

(F = 0) Eq. 3.32 gives a representation of the potential 
cf , anywhere in the volume, in terms of the values of 
CP and 2$- on the surface of the body. The values of 

1 2>^i 

2<(/dn on the surface of the body are given by the boundary 
conditions but the values of cj> are not known. In order 
to solve the problem, it is thus necessary to obtain first 
the values of c j> on the surface. This can be done by 
letting the point P of the volume V approach a point P* 
of the surface. Then Eq. 4.9, with F = 0, yields 


An £ Cf>(P„ , fc ) 

-- £3 ; ' # ^ £ % 1 r t iTT’l 


Al r 


+ 


f [«■ '•(?!- s % 44)]' < 
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In this Appendix, it is shown that, in the limit, Eq. 4.9 
is still valid if the definition of the function E is 
generalized as follows 


E = 1 
E = 1/2 
E = 0 


outside ll 
on A 
inside 


C.2 


By letting P — *P*, the integrands become singular in the 
neighborhood of P* . Thus, it is convenient to separate 
the contribution of a small neighborhood* of P* , which 
will be indicated as . Thus, Eq. 3.32 can be re- 

written as 


4/l£ = - |{ O.s-v.f - 

z T _r 


1 41 


d<£ ] r J_ dZ r 

J r p {V,S T / 



1 4£ 4 /j_\ ] >p T dz T 
<£ dt, dt,\ r p/ J 1 |^S T | 



-i Ai f it 

a; c*jt, \ 



dZ 1 

|7,S r | 


+ C.3 

*The neighborhood Z! t is a small circular surface element, 
with center P* and radius £. 
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where is the contribution of the neighborhood of 

P* , given by 



//b 




i dS 

~dt, dJb. 
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r P JV,S T ( 
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dt> 
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T dlS 
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C. 4 


In order to simplify the discussion of the limit of J t 
as P — * P*, the steady incompressible case is considered 
first. The results are then extended to the unsteady com- 
pressible subsonic case. The supersonic case is now under 
consideration . 

C.2. Steady Incompressible Flow 

For steady incompressible flow, Eq. C.3 and C.4 reduce 


to 


4n£ cj> - 


_ [( _2t 1A1 


1-lr 


-Hum 

Z-?t 


dZ 


C.5 


and 



C.6 


116 



The analysis of cT £ is highly simplified by using local 
coordinate X, Y, Z, with Z normal to the tangent plane 
(directed from E = 0 to E = 1) . Then, separating terms 
of order £ , Eq. C.6 reduces to* 


= -(%\ [[ , ■■= - 

V3Z. L JJ '/x, i rY, s >Z 1 


*, 2 + 1', 2 


- if II 2. — 1 - MX.dV. + 0(e) 

U /^~Y 7+Z 1 




C. 7 


where the subscript * indicates evaluation at Pa- 
using polar coordinate 

R = V * V? 

0 = icm ' V 1 - . 

one obtains 


£. 
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By 


C . 8 
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♦Note that X = Y 


o, z x = 
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Noting that 


j_ iL f 1 \ , .(eN-z 2 ) 2 = i 2 f-— i= ) 

z dz v/i?s z*/ K Wfc z tz*/ c.io 


Equation C.9 becomes 



= -2n 
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- z " t z 


i/r^Tz' j 


+ O(t) 
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Finally, by letting P go to P* , (that is, Z -► 0) , one 
obtains 



2 * 


t 


t 0(e) 


C.12 



where the upper (lower) sign holds for Z > 0 (Z < 0) , that 
is, when P originates outside (inside) the surface TL ; 
correspondingly, the function E assumes the value E = 1 
(E = 0) . 

Finally, using this result in Eq. C.4, one obtains 





9c p 

0 n, 
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Mi 


Ml + 0(i j 


C.13 


I -I. 




Note that, in both cases (P inside or outside l ) , 



C . 14 


Furthermore, r* is the distance between the dummy 
point, P 1 , and the control point (on the surface Z. ), 
P*. Hence, by letting £ go to zero, Eq. C.13 yields 


E . <f, - § 




1 Ml 
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§ ikd) 


Ml 


C. 15 


It should be emphasized that the limit L-+ 0 is now per- 
formed with P on the surface Z . This implies that the 
contribution of Z £ is now of order £ . In order to 

clarify this point, consider the quantity 
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L 

£. 


— f . ' 



and note that 


C.16 


j U I £ 

- Ii'm 

JL. - *_ 

-- Sig 'I (2) 


l «.-0 

1-0 

l 2 l 


C.17 


whereas 


j li«n \ : l'*i 0=0 

U-° £ ‘J «.-° 


C.18 


The difference between these two limits is due to the fact 

. . . T £ 

that, in the limit (as Z — * 0) , the integrand of 1^ 
behaves like a Dirac delta function and hence, its con- 
tribution for a domain which excludes the singular point 
is zero. 

It may be worth noting that, in Eg. C.12, the sequence 
of limits indicated in Eq. C.17 must be performed, whereas 
in Eq. C.15, the one indicated in Eq. C.18 must be used. 
Finally, it is shown that the results obtained here 
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are equivalent to the definition of E given by Eq. C.2: 
note that for steady incompressible flow, Eq. 4.9, with 
F = 0, yields 


4/1 E f = ' I IS; 4 


C . 19 


Note that Eq. C.19 must be used if P is outside or 
inside the surface, whereas Eq. C.15 must be used if P is 
on the surface. However, by comparing Eqs . C.15 and C.19, 
it is easily seen that Eq. 19 is valid everywhere (outside, 
inside and on the surface 2, ) , if the convention is made 
that E is given by Eq. C.2. 

C.3. Unsteady Compressible Subsonic Flow 

In order to simplify the analysis of unsteady sub- 
sonic flow, it is convenient first to analyze Eq. C.4 with 
the nonrestrictive assumption that the frame of reference 
is connected with the undisturbed air; this implies that 


17 = O M*. o 




C . 20 


This is considered in Subsection C.3.1. The general case 
is discussed in Subsection C.3.2. 
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C.3.1. Frame of Reference with U 


0 


With the assumption C. 20,* Eq . C.4 reduces to 
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ll 
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dl 
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C. 21 


where, according to Eqs. 2.38 and C.20, 

T = — C . 22 

In order to examine Eq. C.21, it is convenient to 
introduce a time shift so that t = 0. Furthermore, let P* 

be a point located on the surface at t = 0; then it is 
convenient to consider a frame of reference with origin at 
P* and Z axis directed along the normal to the surface of 
the body (Fig. C.l). 


* 

Note that, strictly speaking, this assumption makes Eq. 1.3 
meaningless. Hence, this section should be considered as 
a mathematical introduction to Subsection C.3.2. Physically 
meaningful results can be obtained if <f> is replaced by <£ 
in this whole Subsection. 
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Thus, t;he equation of the body is given by 

S, z.- T(x 0 v., t.)= o c . 23 

with 

T (0,0,0), o 

1 1 (0.0,0)-- 0 

|i ( 0 , 0, o), o C.24 

but, in general 

( 0 , 0 , 0 )* C t 0 c.25 

2fc, ' 

By using a Taylor series expansion and neglecting higher 
order terms, Eq. C.23 reduces to 

S = Z f - ct ( -- 0 C . 26 

where Eqs. C.24 and C.25 have been used. Eq. C.26 shows 

that c is the velocity of the surface of the body at P = P* 

and t = 0. Hence, in an infinitesimal neighborhood of 

T 

P* , the surface S is given by 

S T = Z,-cl + £ ^ + V,M z,-*T 

C . 27 

= z ,,-ct + X + = o 
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with 


R = \J X,% F C . 28 

Next, it is important to note that, for P = P* (i.e., 

Z = 0 and t = 0) , the surface <2 does not have a tangent 
plane in P* , but rather a tangent cone. In order to see 
this, consider the normal to given by ^7, S T /|V, S T [ 
where, according to Eq. C.27 


5 T k t l v { /F 71* 


C . 29 


where K is a unit vector in direction Z. Equation 

C.29 shows that the surface H T = 0 has the shape of 

« 

^ T 

a cone in an infinitesimal neighborhood, , of P*. 

L w 

For, Eq. C‘. 29 can be rewritten as 

V.S r K 


with 


H = i 


C.30 


C . 31 


In Eq. C.30, K is constant whereas H is directed 
along the tangent to the surface Note that K + H 

is directed along the normal to This configuration 

is sketched in Fig. C.2 where K and H are shown at two 

I »» 

points, P^ and P 1 , in the neighborhood of P* . 

In order to analyze the contribution it is 

ZT 
£ ' 
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defined as the plane normal to K at the point P* . Note, 
however, that for Z = 0, the surface ZL T is. a hyperboloid 
with axis parallel to K. 

Noting that 


dl T dl r 2nRdR. 

|^s T | ~ 1 3s T /az ( | i + i 2 a 


and that the first integral in Eq. C.21, as well as T, 
is of order e. , Eq. C.21 reduces to 




In order to evaluate the integrals, it is convenient to 
use, as variable of integration, the variable 


z, . z, - it 


C.34 


Note that, according to Eg. C.27 


Z, -- - i-r L 7 - z f 

o« 

where 

z = Z - ct 

Note that, according to Eq. C.27 

c R JR 


('*} -tK ■ • i £ 


By using Eqs . C.35 and C.37, Eq. C.33 reduces to* 
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C.35 


C. 36 
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■ n cp(R,,t) 
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+ 0(C) 
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Z,(t) 

2nif (?,,*) ± I £liL dZ, 
«- I z* 


Z,(o) 


* 

Using Eq. C.32, C.35 and C.37, it is easy to show that, 
as anticipated, the first integral in Eq. C.21 is indeed 
of order £ 
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C . 38 
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fcr 0 


where Z^R) is obtained by solving Eq. C.25 with respect 
to Z^. This yields 


l- -^f,[ r2 W2* + 0-r*;^ 


C . 39 


with r = c/a. Hence 
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C . 40 
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and 
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since, according to C.36 

_L dz c _ _ r 

a* d t a o» 

Combining Eqs. C.38, C.40, C.41 and C.42 yields 
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C.41 
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C . 43 


C . 44 


C . 45 
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Finally, combining Eqs . C.43 to C.45 yields (note that, 
for t = 0 , Z ^ = Z ^ and Z = Z) 


= 2/1 cfi?,, o) 


l ! z I /zSti-r 2 )^ 


f ( p - 0 )? 


z f,/ r +z/Jzso-r*n» 

\ r + z/)z| 


C . 46 


and, by letting Z go to zero 

i'* <5‘ t = 2 /t (p (P 4 , 0) ) 


C . 47 


which is identically equal to Eq. C.12. This implies 
that the results obtained in Subsection C.l (in particular 
Eq. C.2), are valid also for unsteady compressible sub- 
sonic flow with U= 0. The case 0 is considered 
next. 

C.3.2. Frame of Reference with tZ. ^ 0 

In Subsection C.3.2, it is shown that the results 
obtained in Subsection C.2 are valid also for unsteady 
compressible subsonic flow with frame of reference 
connected with the undisturbed air ( =0). In this 

subsection, this last restriction is removed. In this 
case, noting that the first integral of Eq. C.4 is of 
order £ , and T is also of order £ , Eq. C.4 can be 
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rewritten as 



In order to analyze Eq. C.48, it is convenient to use the 
Prandtl-Glauert transformation. By using Eqs. 3.37 and 
3.40, Eq. C.20 reduces to 
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C. 50 


Combining Eq. C.49 and C.50 and generalizing the Prandtl- 
Glauert transformation by introducing 


yields 
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C . 52 


+ 0(t) 


This equation can be considerably simplified by noting 
that 


131 



V 5 T = 


v..(s'j -.{z,S-2j.r.,r]' 

V. s - f* (»?**'•)[ 

fe.S- li Mt) . as ^ r J T 

[> 3fc>, J 

- [ 5 - 

|5| (fc. . H.)] t 


C.53 


where i is a unit vector in x Q - direction and 

D = V 0 S - £i X 

*- z>t», 

- $./ |DJ 

«■ > - g- V - r - /i S| = ^ V- 

O to, 7 


with 


^ - 


3S 1 


C.54 


^ w 

Note that Eqs. C.53 and C.54 are in full correspondance 
with Eqs. C.30 and C.31. Combining Eqs. C.52 to C.55 
yields 


C.55 


*‘■*4 




• u { * m 


C.56 
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Equation C.56 is formally equal to Eq. C.33. Hence, 
the results obtained for 0 (Subsection C.3.1), 

are valid for the. general case as well. 

Finally, an important remark should be made. As 
shown by Eq. C. 38, each of the two integrals tends to 
infinity as Z goes to zero. However, their difference 
tends to the finite value given by Eq. C.47. Hence, 
the numerical integration of the two integrals must be 
very careful. 
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Fig. C.l The surface 2 . in the neighborhood of P 
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APPENDIX D 


UNSTEADY WAKE AND LIFTING SURFACE THEORY KERNEL 


D.l. Introduction 

In this Appendix, an explicit expression for the 
evaluation of the unsteady wake (Subsection D.2) is given. 
Finally, for the sake of completeness, an explicit ex- 
pression for the kernel function of the integral equation 
given in Ref. 8 is derived from the results obtained in 
this appendix. 

D.2. Unsteady Subsonic Wake 

In this subsection, a treatment for the wake of an 
oscillating wing in subsonic flow is derived. The same 
assumption made for the steady state wake (Subsection 6.3) 
is made here: the wake is approximated by straight vortex- 

lines, parallel to the direction of the flow, emanating 
from the trailing edge of the wing. For the sake of 
simplicity, the trailing edge is assumed to lie on the 
plane = 0. 

Under this assumption, Eq. 4.19 reduces to 



ft -ttjT 

jj> Q n « ± dz 

z, '> 


+ § $ JL (^J + I* 


D.l 
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with 



In order to take advantage of the fact that there is no 
pressure jump across the wake, it is convenient to follow 
a procedure similar to the one outlined in Subsection 5.2 
By using Eq. 5.11, Eq. D.2 may be rewritten as 

b/2 ao 


'"'O 








j 



-k u 


D. 3 


where 


VI 





dim, 






TS 





D. 4 


137 



where integration by parts has been used and 

'/ 2 


r 

A, 




Z 

X) + 


■ir 


D.5 


Eq. D.4 can be simplified considerably by noting that 

= 0 on the wake: this implies that the last term is 

identically equal to zero. Furthermore, the upper con- 
tribution of the first term is zero because the first 
integral is equal to zero for x.^ = . Finally, note 

that, according to Eq. 5.11 


( 1 UJ 

-_L e Ap( 


“ f cJ x rg 
e U “ 


D.6 


Hence, Eq. D.4 reduces to 



where 


A = 


jc - A, , 


and 


r 

A 


n i/a 




D . 7 


D . 8 
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Note that, according to Eq. 2.38, with x - x^ = -A , 

A. _T = A- - _1 fr . ma] 

<4 ] 



Finally, combining Eqs. D.3, D.7 and D.9 yields 
b/z 

-bit 


D.9 


D. 10 


with (note that <3/02, - - 3/^2 ) 


(x T M ) = - e U - 1 1 

e “r 

1 ft 

L r A 


-6k 


<k\ 


D.ll 


Note that, according to Eq. 5.15, the Kernel function of 
the lifting surface integral equation is given by 


k 


2 


J w 


D . 12 


D.3. An Explicit Expression For J w 

In order to derive an explicit expression for J w , 
consider the integral 
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<k\ 




A 

In order to evaluate consider the classical 


Mr, - A . M JamF - A , 

p>R 

with 

R = jS / 

Note that 

Jl + U* , _L jp 2 £* + M l 0% R z j - 2 M Aj/vTi 1 + A* 

(sR L 

* _L £(fcSA*) -2MA iT AS R* +M 2 -A l J ,/2 ’ 

= feg 2 -MA 

and 


du 

<*A 


_L Fm_a | 

jaR L fo+k 2 - . 


I A Mg 1 - MA _ __ \/7+h«- 
(3 R r A r A 


D. 13 


transformation 


D . 14 


D . 15 


1/2. 


D . 16 


D . 17 
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Hence, using the transformation D.14, Eq. D.13 becomes 


A 

J = 

w 


2 


i cjR 

f u - 


<Am 


D.18 




u, 


where 


M.= 


M /fx-*,) + 2 Z _ Of -x , ) 



Consider the integral in Eq. D.18 


1 


oO 

r 




-ixu 

e — <Am 

'/TT** 


where 


D. 19 


D.20 


k ■ -Si-R » d . 21 

A 

In order to evaluate I it is convenient to make use of 
the contour integration. Consider the integral 



C 


where C = C^ + C 2 + + C 5 is the contour indicated 

in Fig. D.l. The point + i and - i are branch points: 
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by connecting them by the branch line shown in Fig. D.l, 
the integrand can be treated as an analytic function; 
thus, the Cauchy's. Theorem yields that 

I, + + *3 + *4 = 0 D. 23 


where 


{ i/T+w* 


(Am 


D . 24 


On the other hand, if the radius of the circle goes to 
zero and the circle of goes to infinity, one obtains 


0 


D.25 


and, by the Jordan lemma 



Furthermore , 

I, — * 1 = j £ <Am 

«, VTTu* 


D.26 


D . 27 


whereas 



- 1 


D.28 
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and finally (Ref. 23) 


T 5 — h = 


■iXK 


e <ui 


V7T, 


-X^ 

^ AaT 

^*7 


- K © 

* de 


KW 


D.29 


where K 0 (>0 is the modified Bessel function of second kind 
of zero order. 

Hence, by using Eqs . D.25 to D.29, Eq. D.23 yields, 
in the limit 

Ol 

A 

I = 

« 

U, 


- ,KM A A . , -<•'»<« 

£ flu - _ - Xr - (Am 


I + u 2. 


Vl+ U* 


D. 30 


Next, this result can be used in Eq. D.18 to yield 


A A 


-*Xwl 




JL : 2! : - K (X) i3L _ 1 1*. « e — fitw - 


dt 


Ql 02 


V777* 


Vl + M, 1 


0u, 

32 


D.31 


where the relation (Ref. 23) 


3 fc(*) . . K (x) 

a * 


has been used, where is the modified Bessel function of 
second kind of first order, given by (Ref. 23) 
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K (*) -- (t !,(») * -t 

Z _! f(l.i . )-i J-lfej 

2 «=o «,!(»»()/ l 2 ">*' ; 2 

where 7 = .4772157 is the Euler constant. 
Finally, the indefinite integral 


F M = - 1 H 


.iKM 

a - a* 
Y 1+ iF 


D. 33 


D.34 


is analyzed in Subsection D.4. The results are given by 
Eqs. D. 41 to D. 44. By combining Eq. D.ll, D.13 and D.34, 
one obtains 


(f, - *) A 

i (».,*)= - e ““ O'”*) 


CJW 

u 


«Jd- (x,-x) 


_L [ F (-0 - F(u,) 

K dl 1 


. i'*t u. D . 3 5 

- -f l^L - 

V I + w 2 ^2 


where H is given by Eq. D.21 and 


go SR _ ^_|_9R _ A 2 x ill. . i± -±lL 

” _ . *-■% / ^i 3 ' 
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, u, tfz R. se R l U m V(ij-q l )V(e"-? l ) 1 


D . 36 


since, according to Eq . D.15 


1E_ - 

32 


-P 


2 - e. 
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-p 5 


2 - 2 i 


D . 37 
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and,' finally , according to Eqs. D.16, D.19 and D.3 8 




(Z-2, )(*,-*) J_ 

r p 


Note that, as shown in Subsection D.5 (Eq. D.63) 

f i-g * - f- * i>> D . 

where I^(rO is the modified Bessel function of the first 
kind of first order (Ref. 23). 

Finally, combining Eqs. D.35 to D.39 yields 
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('l-1.) 2 +(2-2 1 ) 1 r p 


(H-Hf +(«? -2«) 2 


F(uJ 


D. 40 


(?-*,) 


As mentioned above, an explicit expression for F (u) 
is derived in Subsection D.4, where it is shown that 


F(u] . L f» 


D. 41 


where F (u) can be evaluated by using the recurrent formula 


Fjwlr JL (- 1 )" m.” /|7? u" * + !*.,(«) 

"I"' 2 ) 


• n 


D.42 


with x given by Eq. D.21 and 


^ (u) . i« i/| + M 


and 


T\H = --£• 

2 . 7 


U J I ■+ M l _ (u •+ Jl4 U l ) 


D. 43 
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D.4. The Kernel of the Lifting Surface Theory 

In this Subsection, the results obtained thus far 
are used in order to derive an explicit expression for 
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the Kernel function of the lifting surface theory. In 
order to do this, it is convenient to rewrite Eq. D.40 


as 


where 


At 


AS At 

V -U-z,)K 


D. 45 
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According to Eq. D.12 
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and 
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In order to show that this expression is equal to the 
one given in Ref. 8, it is sufficient to note that 


A A 


F (uj - F(-<) - p* e <U + 


= -tX I Jii dx - j x f i K n f J (*) - L (k) 

1 SJ77 z L • ' ‘ J 

where Eq. D.66 has been used. 

Finally, combining Eq. D.48 and D.49, yields 


iv u v • ' 
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\/ I 4 U* 




/(>,-«)*+ 


in agreement with Eq. D.8 given in Ref. 8. 
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D.5. Evaluation of the Integral F 

Consider the indefinite integral 

• < X M 

M £ dM 

V i + n 1 - 

This integral can be integrated by series as follows: 



A 

F 


M. 


J Am = Z f n 








(*»-•)! 


D.52 


with 



D . 53 


Note that the interchange between integration and summation 
signs is allowed by the uniform convergence of the series 
of the exponential function. 

The integral given by Eq. D.53 can be evaluated by 
using the recurrent formula* 



«• I 

u 



- (*-•) 


D.54 


In particular, disregarding the constant of integration, 



*Note that, for u > 1 , f is divergent as n goes to infinity. 
However, the ratio $/(«-»)! (where f is evaluated analytically) is 
convergent to zero. 


149 



Hence, Eq 


D.56 


D . 52 may be written as 



where 
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or 
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with 


F, s . ; * fiTZ 
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Note, in particular, that F^- i> " 0 and 


Hence, F (-i) = 0 for n = odd and 
n 


C-v)= _iL F U) . -ni(±) L 




for n = even. Thus 
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Note also that F 2 (0) = 0 and 


F (o) = -tH 


Hence F (0) = 0 for n = even and 
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. , >o i 
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(zl r(H*)r(a) 

for n = odd (where T is the Euler gamma function or 
factorial function) . Thus 
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where L ^(H) is the Struve function (Ref. 23). 

Combining Eqs. D.63 and D.66 yields 

- i 

a r 

f (, t )-F(°j = -. K J±s dM T-ni \ l(*i_ L»L;« 

I 2 L J 


in agreement with the well known relation (Ref. 23) 
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branch Is 

line- " " T 



The contour of integration 



APPENDIX E 


SOURCES AND DOUBLETS ON A TRAPEZOIDAL ELEMENT 

E.l. Trapezoidal Element 

In this Appendix, the effect of the sources and 
doublets distributed on a trapezoidal planar element are 
obtained in analytic form. As mentioned in Section 6, 
it is of interest to consider planar elements described 
by the equation 

2, -2, e = *(*, * *ic) -f fOf.-tc) E.l 

where the subscript c stands for center of the element. 

The boundary of the projection of this element on the 
plane z = 0 are given by 

^(Vl) * *, * \ 4 

E.2 

U ' t * H„ 

Equation E.2 represents a trapezoid and the element 
defined by Eqs. E.l and E.2 is called trapezoidal planar 
element (see Fig. E.l). 

E.2. Doublet Distribution 

Consider first the integral of a doublet distribution 
of unit density over a trapezoidal planar element, given 
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by 


• He. (-:) 

t P 

- - / h jl£f'[M ; *£<*•-« 


■i Ac, 


E. 3 


.5 5 
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*(*.-*) + f (4,-f) -(*, - 4 ^7i 
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where 


9$ /gg, 

I as /2a, 


on upper surface 
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E. 4 


and use has been made of the fact that, according to 
Eg • E • 1 

£ 8 - 21 111 . 

2*1 / a«, 

E. 5 

8 / 35 

r ~ / ae, 

It should be noted that, for trapezoidal planar elements, oc 
and p are constant and that, according to Eq. E.l, 

1 = (a, - 2 ) - *(*,-«) -ft ('i.-'i) 

E.6 

--be - *)-«(*■<-*) - } farl) 
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Thus, Eq. E.3 reduces to 

% *, P + d tforl) 

A f . 

1 :.J UJ 

D J J 

U * < (t - i) 


1 M, 


E . 7 


with 
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(*r*T + (Hr 'if * (s + -c *,-*) + 


E.8 


This integral can be simplified considerably if the new 
variables are introduced: 
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These new variables represent coordinates in a new frame 
of reference with origin in the point x, y, z and base 
vectors 


- 1 ) 
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It should be noted that the vector k is normal to the 
plane of the element. It is also immediately verified 
that i, j, k are mutually orthogonal unit vectors and 
thus represent a (right-hand) reference system. 

Note that z.^ is not an independent variable: com- 

bining Eqs . E.l and E.9 yields 


( I -w z )(x ( -x)+ ipfarl) ( l+ * 4 ) 


1 S 


(•+** ♦M(h.-h) + 


(l + *‘P ( NS '*?)' 1 


% -- l 0+* l +p a )' 


E.ll 


where % is a constant given by Eq. E.6. 

Using ^ and i-j as new variables of integration, 

Eq. E.l reduces to (note that dx^dy^ = (1 + "of? + T d^d * 7 ) 
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P| p -. [(uv‘ ’ 

^ [( i ♦-«* 1 ^ I] Ci «*)’’ (' * ' '■’* f*)’ ' 1 


E. 13 


157 



and 
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with 
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e<^ -t- (l ■*■«* ) d. 
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Integrating Eq. E.12 yields 
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Consider the indefinite integral 
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1 J y«i i L J 

This can be integrated by standard methods of integration; 
using the transformation 
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and integrating yields 
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E.22 


and, returning to the original variable, 
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Finally, by using Eg. E.22, Eg. E.18 reduces to 
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with 



•A >\ 
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* §,m *)m 


E.26 


E.3. Source Distribution 


Finally, consider the integral of a source distri- 
bution of density over a trapezoidal planar 

element, given by 


Ji r 2*i 
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Using the transformation introduced in the preceeding 
subsection yields 

\ l,ti> 
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and, by integration 
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Consider the indefinite integral 
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Combining Eqs. E.31 and E.32 yields (note that 4, *\ - 4 - - ) 
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Combining Eqs. E.29, E.33 and E.34 yields 
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e trapezoidal element in the 
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APPENDIX F 


INITIAL CONDITIONS 


In this Appendix, the generalized Huygens' principle 
is further generalized to cover also the case in which 
the problem is formulated with initial conditions. In 
this case, the function E, defined by Eq. 3.1, must be 
assumed to be zero also for t - 0. Thus, a branch of 
the surface is given by the equation* 

5 = t, = 0 F . 1 

For this surface, VS = 0 and 3S/3t ( = +1. Furthermore, 
S T = t - T and |vs T | = | V T j . Thus, Eq. 3.32 can be 
generalized by adding to the right-hand side the term 


4n E f 
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at 
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21 
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IVT | 


* 

The condition that O S is directed like the four- 
dimensional outward normal v is satisfied (see Fig. 2). 
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As shown in Appendix B, for M = 0, Eq. F.2 reduces to 

9 

the classical Poisson formula (Cauchy problem for the 
wave equation) . 

Similar results can be obtained for the supersonic 

flow. 
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APPENDIX G 


SMALL THICKNESS FORMULATION 


G.l Introduction 

As mentioned in Section 6, the integral equation for the 
potential ij) has the disadvantage that it becomes singular as 
the thickness goes to zero. In this Appendix, the charac- 
teristics of this equation are analyzed more in detail. In 
particular, the behaviour at the trailing and the leading 
edges (where one might consider that the thickness is zero) 
are examined (Subsection G.2), vortex layer interpretation 
is given (Subsection G.3) and an alternative approach to 
solving the problem is suggested (Subsection G.4). 

G.2 Leading and Trailing Edges 

Consider first the trailing edge. In order to obtain 
the integral equation at the trailing edge, it is convenient 
to use the same procedure used in Appendix C, that is, to 
obtain the integral equation as the limit (when the control 
point approaches the trailing edge) of Eq . 6.1, or for simplicity, 
the corresponding one for incompressible flow 



For the sake of simplicity, the surface 2 is replaced by the 
smoothed surface 21 in an infinitesimal neighborhood of the 
trailing edge (see Fig. G.l). 
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Next, the limit value of Eq. G.l as the control point P 
goes to P* , is considered. Following the same procedure used 
in Appendix C, one obtains immediately 



where (^ TE |) the u PP er (lower) value of (j> at the trailing 

edge. Similarly, if P approaches P* from the bottom, one 
obtains 


f£,U dZ j G - 3 


h 

However, the two equations are not independent since both are 
equivalent to 


{**£«-* 


G. 4 


Note that in general for lifting bodies (f) ^ f et ’ T ^ us ‘ 

one has only one equation (Eq. G.4) for two unknowns ^ and 
f ’ Hence / an additional condition, the Kutta condition, 
must be given in order to make the problem complete: 


^ C f,T£ = (^>.h ” C?.e ) 7 _ = ( „ difif \ _ q G.5 

!e dit c>Jt he 

The implication of this on the numerical formulation described 
in Section 6 are obvious. If the control points of the upper 
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and lower surface are very close to the trailing edge, the two 
equations are both very close to Eq . G.4 and thus, the deter- 
minant is very close to zero. Hence, use of small boxes in 
the neighborhood of the trailing edge yields elimination of 
significant figures even for thick wings. An alternative 
formulation is then required (Subsection G.3). 

It may be noted that the above problem does not exist 
at the leading edge since there , the upper and lower values 
of the potential are equal. Incidentally, however, it may 
be worth noting that, as is well known, at the leading edge 
the hypothesis of small perturbation fails. For, at the 
leading edge, the boundary conditions are given by 



9S/0X _ y\ y - i 

- i K 

|7S| 


G. 6 


(since the normal n is almost parallel to the x-axis) , which 
is in contradiction to Eqs. 4.5 and 4.6. It may be noted 
however, that the boundary condition is still used in its 
exact form and that the terms neglected in Section 4 are 
still of the same order of the nonlinear term of the differ- 
ential equation. Moreover, this yields only a local effect 
which eventually can be analyzed with the method of matched 
asymptotic expansions (see for instance, Ref. 24 ) . in any 

case, the error involved is smaller than the one obtained in 
the lifting surface theory where the value of /0X is 

infinitely large. 
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G.3 Vortex Layer 


As is well known, a discontinuity in the potential 

corresponds to a vortex layer. The direction of the vortices 

in the layer is given by the direction of the line Af= constant 

3 

and the intensity of the vortex is given by variation of • 0S ~ 
where s is the direction (in the layer) normal to the "constant 
lines", that is, the direction of the perturbation 
velocity component in the plane of the vortex layer. Phy- 
sically speaking, the doublet integral can be interpreted as 
a (zero-thickness) boundary layer (note that, with this inter- 
pretation, E = 0 inside Z should be replaced by C^= 0 inside 
Zj , which implies that the perturbation velocity is identi- 
cally equal to zero inside Zj ). Finally, it may be noted 
that this implies that this formulation does not yield any 
phenomena of the type encountered in lifting-surface formu- 
lations (with elements inclined to the flow) , for which "wakes 
emanating from points near the body leading edge will thread 
through the body surface near its trailing edge" (Ref. 25 
p. 446 ) . It should be noted that the tangent plane approxi- 
mation (Section 7) yields a similar phenomenon. However, the 
tangent plane approximation should be considered only as a 
numerical approximate procedure (with controlled error) to 
solve a physically well posed formulation. 

G.4 Alternative Formulations 

Finally, it is worth noting that alternative numerical 
formulations can be used in connection with the theoretical 
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formulation presented here. Again, for the sake of simplicity, 
the discussion is carried out for the steady incompressible 
case. Consider first the normal derivative of the basic equa- 
tion 



By taking the limit as the control point approaches Z_ , one 
obtains an integral equation (different from Eq . 6.2), which has 
the advantage that if the thickness is equal to zero, the opera- 
tor is not singular: the limit equation is the one used by 

Haviland . ^ 

Two other alternative formulations are obtained by intro- 
ducing a convenient flow field inside the surface ZL (see 
27 

Lamb ) . If the value of the potential is continuous across 
the surface, one obtains 



G . 8 


where S is the intensity of the source distribution (equal to 
the discontinuity in normal velocity). On the other hand, 
if the normal velocity is continuous across the surface, then 
one obtains 



G. 9 


where D is the intensity of the doublet distribution (equal 
to the discontinuity of <J> ) . By taking the normal derivative 
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of Eqs • G.8 and G.9 and imposing the boundary conditions, 
one obtains two integral equations for the unknowns S and D, 

respectively. These methods are used, for instance, by Hess 

28 6 
and Smith and D jo jodihard jo and Widnall, respectively. 

It may be noted that more complex formulations can be ob- 
tained by combining two or more of these four basic methods, 
Eqs. 6.2, G.7, G.8 and G.9. The advantages and disadvan- 
tages of the four basic methods are briefly discussed here. ' 
Note first that the "source method", Eq. G.8, is limited to 

nonlifting bodies: extensions to lifting configurations must 

2 9 

include doublet distributions as well. On the other hand, 
the other two methods, Eqs. G.7 and G.9, if applied to a 
closed surface, involve a singular operator (which yields a 
determinant equal to zero). For, according to Eq. 6.6, one 
obtains, outside Z , 



and thus its normal derivative is zero: hence, Eqs. G.7 and 

G.9 have a nontrivial (constant) solution for the homogeneous 
problem.* Thus, the method presented here is the only one of 
the four basic methods, which can be applied to "closed-surface" 

* 

It should be noted however, that these two methods can be 

2 6 

used for open surfaces (lifting surfaces ) , or for the analysis 
of the transient response.^ 
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description of lifting conf igurations . On the other hand, as 
mentioned above, combinations of two or more of the four basic 
methods can be usefully employed. Various combinations are 
now being explored. 
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APPENDIX H 


SUPERSONIC DOUBLET 


H.l Introduction 

It should be noted that, as mentioned in Subsection 3.5, 
the supersonic Green function (Eq. 3.44) has an infinitely 
large discontinuity at the Mach cone. Hence, the normal de- 
rivative (supersonic doublet) has meaning only in 

the theory of distributions . ^ 

The scope of this Appendix is to obtain the correct 
definition of the supersonic doublet within the theory of 

functions (as distinct from distributions) . In order to 

22 

obtain this, the Lighthill definition of distribution, or 
generalized functions (as limit of regular functions) will be 
employed. For simplicity, the steady case is considered. 

The potential steady supersonic flow around a body of arbitrary 
shape is described by the equation 


# 




H.l 


(q dl. - dL 

•r T 0KJ ' 

where Cf> is the perturbation potential, Z is a surface surroun- 
ding the body and the wake and (see Eq. 3.44) 


2 ni 


~i 

z 


O 




H.2 


x-x, < 


is the steady supersonic source or Green's function. In Eq. 
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since 


B.2, it is assumed for simplicity that B = ) = I , 

by using the supersonic Prandtl-Glauert transformation, ex- 
pressions formally equal to Eqs . H.l and H.2 can be obtained. 

For the sake of simplicity, an element of planar surface 
ZJ. e , parallel to the x-axis is considered. By suitable 
change of coordinates, this element can be reduced to be 
lying in the plane z^ = 0 . In this case 




0 


0 

Pa 


H. 3 


Furthermore, the origin will be assumed to be moved to the 
point (x,y,0), so that Eq. H.2 reduces to 

<3 , - 


o 




H . 4 


In conclusion, the supersonic doublet integral reduces to 


Zo 


H . 5 


with G given by Eq. H.4. This integral is analyzed in the 
following . 

H.2 Modified Supersonic Doublet Integral 

As mentioned before, following the Lighthill approach, 
the Eq. H.4 will be considered as the limit of a more manage- 
able function. A suitable choice is (see Fig. H.l) 


G, - - - — 


r * ** 


H . 6 


2-n 


r B* £ 
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Furthermore , 

d-O *, H. 7 

Hence, the function has no discontinuity for all the values 
of x, and y. In conclusion, the modified supersonic source G 


is given by 


G ^ * 

Z71 


li 


^ **)’ x *,s - V +2i ) 

*, < CX5 


= O 


Consider the modified supersonic doublet integral 

h - ^ if ^ **• 


H. 8 


H. 9 


where Z, , Z z and Z 3 are shown in Fig. H.2. For the sake 
of simplicity, two edges of 2 0 have been assumed to be 
parallel to the x-axis. Hence, 

A * 

( 

2n 


A / a A % 


H. 10 


where 




and 




Z, 3z * ’*■ 

v 

-- 2 Z J H / 

i // S 


H. 11 


J- d«, 




-- -&hj 




H. 12 




i 

r 5 
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Using the transformation 


X 


i 



42 1 


one obtains 


and 


My, = 


- Si- Mix', 

*t 


H. 13 


H. 14 



L 


y 

. I 4 1 * f _/g^ K, t 4 ~g*_ 

£ 0 


I t4 


+ jr, W ' ] V Y 

1 "l/eV^Vo 


I + 2* uf 'r^JP^WP ) 4 2J_uY. .i- V _ \ 

£ ' £*+4* / I 1 1 m z l'/eVY +z i / 


H. 15 


where 


1 - 2z 


Y / 

f ^ - _ 2 2 I - X — 

*f,‘ u l ^ J l 1 * 1 ' 


H . 16 


is the contribution from the edge x = X(y) 
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On the other hand, combining Eqs. H.12 to H.14 yields 


A 

1 


V 




J 

o 


H j . 


1 -- - ?J_ 

r B e L 




y 

j (t 4 / ) - 


^ cIh, 
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Bv expanding in power of i one obtains 



l + 2i 





4 



I*! 




- 1 


£ V 


ii| /yvT- 


0-T 


r\f 




= I + o(«) 


Finally, combining Eq. H.10 and H.19, one obtains 

T - - _L I t 0(e) 
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H.19 


H . 2 0 


H.3 The Supersonic Doublet Integral 

By letting € go to zero, one obtains 




-A 



H. 21 


where I, given by Eq. H.15, is the contribution of the boundary 
x = X(y^)., Hence, the rule is that the contribution of the 
boundary r^ = 0 (which is infinity if the theory of function 

« • 1 " D 
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is applied carelessly) is zero. 


Note in particular that if X = const. 

2 . L r ; Cz) — y — 

t> 2 n. tin j i l‘/x 2 ~y 2 -^i t H.2 2 

Finally, note that 

r p£ : °) = V X = 3 H. 23 

which is the same result used in subsonic flow to prove that 
E = 1/2 on the surface . 
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